stanEx2
execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}ex5 regression
ex5.stan
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> s;
}
model{
vector[N] m=X*b;
y~normal(m,s);
}
generated quantities{
vector[N] y1;
vector[N] m1=X*b;
for(i in 1:N){
y1[i]=normal_rng(m1[i],s);
}
}normal linear models
n=30
mdl=cmdstan_model('./ex5.stan')single regression
tb=tibble(x=runif(n,0,9),
y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)
qplot(data=tb,x,y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -105.796
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -12.8956 0.000206144 3.95095e-05 1 1 19
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.mle## variable estimate
## lp__ -12.90
## b[1] -0.20
## b[2] 0.99
## s 0.93
## y1[1] -0.83
## y1[2] 8.06
## y1[3] 1.77
## y1[4] -0.74
## y1[5] 7.03
## y1[6] 7.43
## y1[7] 4.32
## y1[8] 6.94
## y1[9] 6.45
## y1[10] 3.48
## y1[11] 8.73
## y1[12] 6.94
## y1[13] 0.73
## y1[14] 7.45
## y1[15] 5.16
## y1[16] 3.37
## y1[17] -0.48
## y1[18] 0.65
## y1[19] 0.76
## y1[20] 6.17
## y1[21] 6.42
## y1[22] 1.03
## y1[23] 8.17
## y1[24] 5.51
## y1[25] 8.03
## y1[26] 3.09
## y1[27] 7.77
## y1[28] 0.37
## y1[29] 8.27
## y1[30] 2.45
## m1[1] 0.39
## m1[2] 6.91
## m1[3] 3.32
## m1[4] 0.28
## m1[5] 6.06
## m1[6] 7.26
## m1[7] 6.42
## m1[8] 5.90
## m1[9] 5.79
## m1[10] 4.45
## m1[11] 7.98
## m1[12] 6.72
## m1[13] 1.55
## m1[14] 7.37
## m1[15] 6.36
## m1[16] 3.20
## m1[17] 0.00
## m1[18] 0.81
## m1[19] 2.79
## m1[20] 5.89
## m1[21] 7.43
## m1[22] 1.11
## m1[23] 8.03
## m1[24] 6.23
## m1[25] 6.85
## m1[26] 1.86
## m1[27] 7.50
## m1[28] 1.12
## m1[29] 8.31
## m1[30] 1.70mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.54 -14.21 1.26 1.05 -17.07 -13.13 1.00 718 868
## b[1] -0.23 -0.23 0.37 0.38 -0.83 0.35 1.01 391 386
## b[2] 1.00 1.00 0.07 0.07 0.89 1.11 1.00 539 749
## s 1.01 0.99 0.14 0.13 0.80 1.26 1.00 1110 772
## y1[1] 0.38 0.40 1.05 1.01 -1.30 2.11 1.00 1188 1852
## y1[2] 6.88 6.88 1.07 1.03 5.13 8.66 1.00 1677 1742
## y1[3] 3.33 3.31 1.02 1.01 1.64 5.02 1.00 1641 1745
## y1[4] 0.21 0.20 1.07 1.05 -1.58 1.96 1.00 1662 1809
## y1[5] 6.04 6.04 1.02 1.00 4.31 7.70 1.00 2009 1801
## y1[6] 7.26 7.23 1.03 1.04 5.67 8.91 1.00 1994 2058
## y1[7] 6.42 6.44 1.06 1.04 4.60 8.13 1.00 1718 1860
## y1[8] 5.92 5.92 1.03 1.01 4.23 7.59 1.00 2131 2009
## y1[9] 5.78 5.78 0.99 0.94 4.19 7.43 1.00 2115 2015
## y1[10] 4.44 4.45 1.07 1.07 2.67 6.14 1.00 2023 1972
## y1[11] 7.99 7.98 1.03 0.98 6.28 9.69 1.00 2061 1968
## y1[12] 6.71 6.71 1.05 1.02 4.98 8.39 1.00 1853 1853
## y1[13] 1.52 1.55 1.06 1.05 -0.26 3.23 1.00 1906 1932
## y1[14] 7.39 7.36 1.09 1.05 5.59 9.18 1.00 2039 1684
## y1[15] 6.36 6.37 1.03 0.98 4.64 8.04 1.00 2102 1927
## y1[16] 3.18 3.16 1.05 1.03 1.45 4.91 1.00 1814 1734
## y1[17] 0.03 0.05 1.09 1.09 -1.77 1.83 1.00 1573 1888
## y1[18] 0.82 0.88 1.12 1.08 -1.08 2.61 1.00 1208 1622
## y1[19] 2.78 2.76 1.04 0.99 1.14 4.55 1.00 2016 1950
## y1[20] 5.90 5.89 1.05 1.01 4.18 7.64 1.00 1991 1856
## y1[21] 7.47 7.45 1.09 1.05 5.67 9.21 1.00 1949 1881
## y1[22] 1.12 1.12 1.07 1.05 -0.58 2.89 1.00 1722 1820
## y1[23] 8.06 8.05 1.07 1.02 6.28 9.78 1.00 2072 1758
## y1[24] 6.24 6.26 1.04 1.02 4.53 7.90 1.00 1963 1997
## y1[25] 6.82 6.82 1.05 1.06 5.11 8.57 1.00 1987 1748
## y1[26] 1.80 1.80 1.03 1.00 0.12 3.44 1.00 1819 1827
## y1[27] 7.49 7.51 1.06 1.04 5.71 9.20 1.00 2016 2010
## y1[28] 1.10 1.07 1.07 1.09 -0.57 2.89 1.00 1583 1806
## y1[29] 8.32 8.31 1.07 1.06 6.63 10.09 1.00 2110 2014
## y1[30] 1.71 1.70 1.03 1.02 0.07 3.39 1.00 1738 1968
## m1[1] 0.37 0.37 0.34 0.34 -0.18 0.91 1.01 389 394
## m1[2] 6.92 6.92 0.24 0.24 6.52 7.32 1.00 2296 1605
## m1[3] 3.31 3.31 0.20 0.20 2.98 3.65 1.00 508 985
## m1[4] 0.25 0.26 0.34 0.35 -0.30 0.80 1.01 389 394
## m1[5] 6.06 6.06 0.21 0.21 5.71 6.41 1.00 2643 1606
## m1[6] 7.27 7.26 0.26 0.25 6.84 7.70 1.00 2083 1647
## m1[7] 6.42 6.42 0.22 0.22 6.05 6.79 1.00 2611 1581
## m1[8] 5.91 5.91 0.21 0.20 5.57 6.25 1.00 2566 1606
## m1[9] 5.79 5.79 0.20 0.20 5.46 6.13 1.00 2478 1665
## m1[10] 4.44 4.44 0.19 0.18 4.13 4.74 1.00 878 1480
## m1[11] 8.00 7.99 0.30 0.29 7.51 8.49 1.00 1670 1591
## m1[12] 6.72 6.72 0.24 0.23 6.33 7.11 1.00 2433 1616
## m1[13] 1.53 1.53 0.28 0.28 1.08 1.97 1.01 396 641
## m1[14] 7.39 7.38 0.27 0.26 6.95 7.83 1.00 2006 1619
## m1[15] 6.36 6.36 0.22 0.22 5.99 6.73 1.00 2634 1581
## m1[16] 3.19 3.19 0.21 0.21 2.85 3.53 1.00 491 932
## m1[17] -0.03 -0.02 0.36 0.36 -0.61 0.54 1.01 390 386
## m1[18] 0.78 0.78 0.31 0.32 0.27 1.29 1.01 389 593
## m1[19] 2.77 2.77 0.22 0.22 2.41 3.13 1.01 449 773
## m1[20] 5.90 5.90 0.21 0.20 5.56 6.24 1.00 2562 1609
## m1[21] 7.44 7.43 0.27 0.26 7.00 7.88 1.00 1973 1619
## m1[22] 1.09 1.09 0.30 0.30 0.61 1.57 1.01 390 604
## m1[23] 8.04 8.04 0.30 0.29 7.55 8.54 1.00 1651 1554
## m1[24] 6.24 6.23 0.22 0.21 5.87 6.59 1.00 2665 1607
## m1[25] 6.86 6.85 0.24 0.23 6.45 7.25 1.00 2333 1616
## m1[26] 1.84 1.84 0.26 0.26 1.42 2.26 1.01 402 683
## m1[27] 7.51 7.51 0.27 0.26 7.07 7.96 1.00 1929 1603
## m1[28] 1.10 1.10 0.30 0.30 0.62 1.58 1.01 390 604
## m1[29] 8.32 8.32 0.31 0.30 7.80 8.84 1.00 1531 1487
## m1[30] 1.68 1.68 0.27 0.27 1.24 2.11 1.01 398 673
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
multiple regression
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -15999
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 44 -13.1144 0.000133986 0.000327327 1 1 86
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.11
## b[1] -0.16
## b[2] 1.04
## b[3] -0.95
## s 0.94
## y1[1] 0.64
## y1[2] 2.63
## y1[3] -3.29
## y1[4] 2.85
## y1[5] -3.74
## y1[6] 4.86
## y1[7] -2.92
## y1[8] -2.71
## y1[9] 8.13
## y1[10] 1.65
## y1[11] -1.58
## y1[12] 1.49
## y1[13] -4.53
## y1[14] -0.26
## y1[15] -1.97
## y1[16] 2.93
## y1[17] -0.68
## y1[18] 2.83
## y1[19] 1.95
## y1[20] -1.32
## y1[21] -0.38
## y1[22] -5.81
## y1[23] -1.32
## y1[24] -2.67
## y1[25] 7.60
## y1[26] 3.43
## y1[27] 2.96
## y1[28] 5.53
## y1[29] -4.26
## y1[30] 5.03
## m1[1] 1.10
## m1[2] 1.35
## m1[3] -3.59
## m1[4] 3.95
## m1[5] -3.21
## m1[6] 4.05
## m1[7] -2.41
## m1[8] -3.91
## m1[9] 7.81
## m1[10] 1.16
## m1[11] -0.61
## m1[12] 3.02
## m1[13] -4.30
## m1[14] -0.16
## m1[15] -3.05
## m1[16] 2.30
## m1[17] -0.76
## m1[18] 2.14
## m1[19] 2.73
## m1[20] -2.01
## m1[21] -1.38
## m1[22] -4.31
## m1[23] -1.23
## m1[24] -4.30
## m1[25] 5.22
## m1[26] 3.06
## m1[27] 2.51
## m1[28] 4.29
## m1[29] -3.01
## m1[30] 4.10mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.38 -15.04 1.52 1.32 -18.45 -13.61 1.00 623 644
## b[1] -0.10 -0.13 0.54 0.53 -0.98 0.84 1.01 330 317
## b[2] 1.03 1.04 0.09 0.09 0.88 1.18 1.00 547 555
## b[3] -0.96 -0.96 0.08 0.07 -1.08 -0.83 1.01 674 710
## s 1.04 1.02 0.16 0.15 0.82 1.32 1.00 1072 774
## y1[1] 1.09 1.08 1.13 1.12 -0.78 2.93 1.00 1109 1268
## y1[2] 1.33 1.30 1.11 1.07 -0.55 3.18 1.00 1984 1784
## y1[3] -3.52 -3.52 1.09 1.06 -5.31 -1.76 1.00 1941 1835
## y1[4] 3.91 3.91 1.12 1.10 2.07 5.74 1.00 1615 2038
## y1[5] -3.19 -3.18 1.08 1.03 -4.98 -1.40 1.00 1928 1668
## y1[6] 4.02 4.05 1.09 1.06 2.20 5.73 1.00 1728 1656
## y1[7] -2.35 -2.36 1.09 1.05 -4.13 -0.49 1.00 1677 1548
## y1[8] -3.91 -3.93 1.13 1.10 -5.74 -2.01 1.00 1975 1994
## y1[9] 7.77 7.77 1.18 1.12 5.87 9.73 1.00 2031 1853
## y1[10] 1.16 1.17 1.07 1.05 -0.64 2.89 1.00 1885 1961
## y1[11] -0.62 -0.59 1.08 1.01 -2.41 1.13 1.00 1878 1838
## y1[12] 3.06 3.04 1.13 1.10 1.22 4.91 1.00 1869 1889
## y1[13] -4.24 -4.27 1.11 1.05 -6.03 -2.40 1.00 2000 2015
## y1[14] -0.12 -0.10 1.07 1.04 -1.84 1.65 1.00 1769 1822
## y1[15] -3.04 -3.05 1.12 1.12 -4.85 -1.27 1.00 1748 1754
## y1[16] 2.27 2.26 1.14 1.09 0.40 4.10 1.00 1818 1637
## y1[17] -0.77 -0.75 1.07 1.09 -2.51 0.94 1.00 1934 1487
## y1[18] 2.15 2.12 1.05 1.04 0.42 3.86 1.00 1940 1930
## y1[19] 2.68 2.66 1.09 1.07 0.95 4.49 1.00 1871 1665
## y1[20] -2.00 -1.96 1.11 1.08 -3.85 -0.23 1.00 1765 1767
## y1[21] -1.38 -1.40 1.07 1.06 -3.15 0.34 1.00 1925 1770
## y1[22] -4.30 -4.30 1.11 1.11 -6.07 -2.43 1.00 2100 1932
## y1[23] -1.13 -1.18 1.17 1.15 -3.04 0.84 1.00 1414 1401
## y1[24] -4.30 -4.33 1.11 1.09 -6.05 -2.40 1.00 1508 1915
## y1[25] 5.23 5.21 1.13 1.08 3.37 7.10 1.00 1760 1811
## y1[26] 3.01 3.02 1.09 1.05 1.22 4.75 1.00 1957 1791
## y1[27] 2.52 2.48 1.10 1.09 0.77 4.39 1.00 1934 1830
## y1[28] 4.31 4.29 1.14 1.13 2.46 6.16 1.00 1909 1760
## y1[29] -3.01 -3.01 1.10 1.05 -4.84 -1.25 1.00 2071 1973
## y1[30] 4.12 4.12 1.14 1.08 2.26 5.94 1.00 2159 1921
## m1[1] 1.13 1.13 0.42 0.40 0.45 1.85 1.01 343 360
## m1[2] 1.37 1.36 0.26 0.25 0.96 1.79 1.00 485 610
## m1[3] -3.56 -3.56 0.31 0.29 -4.06 -3.05 1.00 947 1374
## m1[4] 3.92 3.93 0.44 0.43 3.16 4.62 1.00 734 836
## m1[5] -3.19 -3.19 0.29 0.27 -3.65 -2.71 1.00 957 1414
## m1[6] 4.05 4.05 0.30 0.29 3.57 4.56 1.00 1766 1382
## m1[7] -2.37 -2.38 0.35 0.34 -2.94 -1.79 1.01 427 454
## m1[8] -3.91 -3.91 0.37 0.37 -4.51 -3.30 1.00 1240 1227
## m1[9] 7.80 7.78 0.48 0.47 7.01 8.60 1.00 1818 1254
## m1[10] 1.16 1.16 0.20 0.20 0.84 1.50 1.00 1963 1537
## m1[11] -0.61 -0.61 0.29 0.30 -1.11 -0.14 1.00 948 1303
## m1[12] 3.03 3.03 0.25 0.24 2.62 3.46 1.00 1571 1248
## m1[13] -4.27 -4.27 0.37 0.36 -4.87 -3.64 1.00 695 791
## m1[14] -0.17 -0.17 0.27 0.27 -0.64 0.26 1.00 1069 1340
## m1[15] -3.06 -3.06 0.41 0.40 -3.72 -2.40 1.00 821 1026
## m1[16] 2.27 2.28 0.38 0.37 1.62 2.86 1.00 711 787
## m1[17] -0.76 -0.76 0.22 0.22 -1.11 -0.39 1.00 2137 1576
## m1[18] 2.15 2.14 0.23 0.23 1.79 2.54 1.00 947 1086
## m1[19] 2.71 2.71 0.33 0.33 2.16 3.26 1.00 988 1100
## m1[20] -1.97 -1.98 0.39 0.38 -2.60 -1.31 1.01 379 356
## m1[21] -1.36 -1.36 0.24 0.24 -1.74 -0.96 1.00 615 835
## m1[22] -4.29 -4.30 0.33 0.33 -4.83 -3.75 1.00 1755 1611
## m1[23] -1.19 -1.20 0.46 0.45 -1.93 -0.41 1.01 346 332
## m1[24] -4.29 -4.30 0.36 0.37 -4.87 -3.69 1.00 1511 1290
## m1[25] 5.22 5.21 0.37 0.36 4.64 5.83 1.00 1027 996
## m1[26] 3.04 3.04 0.30 0.30 2.55 3.55 1.00 1579 1560
## m1[27] 2.52 2.52 0.29 0.28 2.05 3.00 1.00 522 725
## m1[28] 4.29 4.29 0.32 0.31 3.79 4.82 1.00 1123 1085
## m1[29] -3.01 -3.02 0.34 0.34 -3.55 -2.45 1.00 1130 1130
## m1[30] 4.10 4.09 0.30 0.29 3.61 4.61 1.00 1181 1205
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANOVA
tb=tibble(c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)
qplot(data=tb,c,y,geom='boxplot')
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -53.587
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -9.48161 0.00018107 0.000101393 1 1 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.48
## b[1] 0.41
## b[2] 1.02
## b[3] -2.56
## s 0.83
## y1[1] -0.72
## y1[2] -0.37
## y1[3] 1.08
## y1[4] -3.57
## y1[5] -2.43
## y1[6] 0.66
## y1[7] 0.70
## y1[8] -2.74
## y1[9] 1.29
## y1[10] 0.22
## y1[11] 0.52
## y1[12] -1.69
## y1[13] 0.10
## y1[14] -3.04
## y1[15] 0.34
## y1[16] 0.81
## y1[17] 0.63
## y1[18] -0.22
## y1[19] -0.24
## y1[20] 0.78
## y1[21] 0.59
## y1[22] 1.58
## y1[23] -0.75
## y1[24] -0.19
## y1[25] -1.82
## y1[26] 0.39
## y1[27] 0.53
## y1[28] 1.48
## y1[29] -0.08
## y1[30] 0.84
## m1[1] 0.41
## m1[2] 0.41
## m1[3] 0.41
## m1[4] -2.15
## m1[5] -2.15
## m1[6] 1.43
## m1[7] 1.43
## m1[8] -2.15
## m1[9] 1.43
## m1[10] 0.41
## m1[11] 0.41
## m1[12] 0.41
## m1[13] 0.41
## m1[14] -2.15
## m1[15] 1.43
## m1[16] 1.43
## m1[17] 1.43
## m1[18] 0.41
## m1[19] 0.41
## m1[20] 1.43
## m1[21] 0.41
## m1[22] 0.41
## m1[23] 0.41
## m1[24] -2.15
## m1[25] -2.15
## m1[26] 0.41
## m1[27] 0.41
## m1[28] 1.43
## m1[29] 0.41
## m1[30] 1.43mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.80 -11.50 1.49 1.27 -14.65 -10.07 1.00 829 937
## b[1] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## b[2] 1.00 1.00 0.39 0.40 0.39 1.65 1.00 901 1041
## b[3] -2.54 -2.55 0.46 0.46 -3.26 -1.78 1.01 780 994
## s 0.92 0.90 0.13 0.13 0.73 1.16 1.00 1829 1485
## y1[1] 0.40 0.38 0.98 0.95 -1.20 1.98 1.00 1819 1893
## y1[2] 0.41 0.43 0.98 0.97 -1.18 1.97 1.00 2033 1810
## y1[3] 0.44 0.44 0.97 0.94 -1.18 2.03 1.00 1795 2090
## y1[4] -2.10 -2.10 0.99 0.97 -3.68 -0.43 1.00 2123 1929
## y1[5] -2.14 -2.14 1.01 1.00 -3.76 -0.49 1.00 1806 1853
## y1[6] 1.42 1.38 0.97 0.95 -0.14 3.10 1.00 1820 1722
## y1[7] 1.42 1.42 0.97 0.94 -0.23 3.02 1.01 1705 1632
## y1[8] -2.07 -2.08 1.01 0.96 -3.71 -0.40 1.00 1860 1971
## y1[9] 1.39 1.36 0.97 0.94 -0.21 2.99 1.00 1878 2015
## y1[10] 0.42 0.44 0.93 0.90 -1.14 1.92 1.00 1898 1971
## y1[11] 0.44 0.42 0.94 0.92 -1.10 1.99 1.00 1948 2001
## y1[12] 0.38 0.40 0.95 0.96 -1.16 1.90 1.00 1553 1959
## y1[13] 0.42 0.40 0.97 0.95 -1.12 2.04 1.00 1992 1943
## y1[14] -2.13 -2.13 0.98 0.96 -3.68 -0.54 1.00 1716 1868
## y1[15] 1.41 1.41 0.98 0.96 -0.17 2.99 1.01 1984 1824
## y1[16] 1.45 1.48 0.96 0.88 -0.18 3.01 1.00 2149 1973
## y1[17] 1.40 1.42 0.99 0.98 -0.17 3.04 1.00 1961 2015
## y1[18] 0.42 0.41 0.95 0.93 -1.14 1.99 1.00 1872 1792
## y1[19] 0.43 0.40 0.95 0.91 -1.09 2.03 1.00 1885 1748
## y1[20] 1.40 1.40 0.98 0.90 -0.22 3.01 1.00 1920 1885
## y1[21] 0.40 0.42 0.97 0.94 -1.20 2.00 1.00 1710 1784
## y1[22] 0.42 0.40 0.96 0.92 -1.17 2.00 1.00 1980 1891
## y1[23] 0.38 0.38 0.98 0.98 -1.24 1.98 1.00 1892 1740
## y1[24] -2.11 -2.14 1.01 0.98 -3.79 -0.43 1.00 1887 1991
## y1[25] -2.15 -2.16 1.01 1.01 -3.78 -0.44 1.00 1672 1842
## y1[26] 0.42 0.43 0.94 0.95 -1.12 1.95 1.00 1821 1757
## y1[27] 0.41 0.40 0.97 0.93 -1.14 2.00 1.00 1807 1915
## y1[28] 1.41 1.39 0.98 0.96 -0.20 2.95 1.00 1575 1934
## y1[29] 0.38 0.37 0.96 0.94 -1.17 1.94 1.00 1959 1912
## y1[30] 1.40 1.41 0.99 0.91 -0.24 3.00 1.00 1917 1869
## m1[1] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[2] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[3] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[4] -2.13 -2.15 0.38 0.37 -2.74 -1.50 1.00 1103 1304
## m1[5] -2.13 -2.15 0.38 0.37 -2.74 -1.50 1.00 1103 1304
## m1[6] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[7] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[8] -2.13 -2.15 0.38 0.37 -2.74 -1.50 1.00 1103 1304
## m1[9] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[10] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[11] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[12] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[13] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[14] -2.13 -2.15 0.38 0.37 -2.74 -1.50 1.00 1103 1304
## m1[15] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[16] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[17] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[18] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[19] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[20] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[21] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[22] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[23] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[24] -2.13 -2.15 0.38 0.37 -2.74 -1.50 1.00 1103 1304
## m1[25] -2.13 -2.15 0.38 0.37 -2.74 -1.50 1.00 1103 1304
## m1[26] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[27] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[28] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
## m1[29] 0.41 0.42 0.24 0.23 0.01 0.80 1.00 924 1376
## m1[30] 1.42 1.42 0.31 0.31 0.92 1.91 1.00 1442 1591
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANCOVA
tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))
f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)
lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])
qplot(data=tb,x,y,shape=c,size=I(2))
plot(tb$x,tb$y,col=1+lv[factor(tb$c)])
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -264.755
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 42 -10.8678 6.35115e-05 0.000874227 0.8883 0.8883 58
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -10.87
## b[1] 0.41
## b[2] 0.97
## b[3] 1.95
## b[4] -1.61
## s 0.87
## y1[1] 7.29
## y1[2] 5.29
## y1[3] 5.51
## y1[4] 7.13
## y1[5] 4.15
## y1[6] 6.09
## y1[7] -0.12
## y1[8] 5.09
## y1[9] 5.60
## y1[10] 2.20
## y1[11] 9.16
## y1[12] 7.63
## y1[13] 3.24
## y1[14] 9.70
## y1[15] 4.27
## y1[16] 5.05
## y1[17] 0.78
## y1[18] 7.85
## y1[19] 6.65
## y1[20] 10.14
## y1[21] 8.69
## y1[22] 0.07
## y1[23] 8.23
## y1[24] 3.84
## y1[25] 4.40
## y1[26] 9.11
## y1[27] 10.32
## y1[28] 3.01
## y1[29] 11.96
## y1[30] 4.94
## m1[1] 5.16
## m1[2] 4.68
## m1[3] 7.18
## m1[4] 6.21
## m1[5] 3.16
## m1[6] 5.42
## m1[7] 0.90
## m1[8] 5.45
## m1[9] 6.15
## m1[10] 2.63
## m1[11] 8.86
## m1[12] 8.13
## m1[13] 3.83
## m1[14] 9.75
## m1[15] 3.94
## m1[16] 5.96
## m1[17] 0.55
## m1[18] 7.32
## m1[19] 7.38
## m1[20] 9.94
## m1[21] 9.64
## m1[22] 0.26
## m1[23] 7.39
## m1[24] 3.09
## m1[25] 4.46
## m1[26] 9.23
## m1[27] 10.20
## m1[28] 1.08
## m1[29] 10.42
## m1[30] 5.88mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.71 -13.37 1.70 1.50 -16.92 -11.58 1.00 752 1036
## b[1] 0.45 0.46 0.41 0.40 -0.23 1.11 1.00 863 945
## b[2] 0.97 0.97 0.08 0.08 0.84 1.09 1.00 1275 1273
## b[3] 1.95 1.94 0.43 0.43 1.24 2.66 1.00 854 929
## b[4] -1.60 -1.64 0.58 0.57 -2.50 -0.61 1.00 640 836
## s 0.98 0.97 0.14 0.13 0.78 1.24 1.00 1526 1177
## y1[1] 5.16 5.14 1.10 1.09 3.31 6.98 1.00 1901 1868
## y1[2] 4.67 4.68 1.03 1.00 2.99 6.31 1.00 1739 1895
## y1[3] 7.13 7.13 1.07 1.03 5.38 8.93 1.00 1814 1936
## y1[4] 6.17 6.19 1.12 1.07 4.32 7.98 1.00 1788 1848
## y1[5] 3.14 3.13 1.02 0.97 1.51 4.83 1.00 1685 1830
## y1[6] 5.39 5.41 1.06 1.03 3.62 7.07 1.00 1704 1740
## y1[7] 0.93 0.94 1.05 0.99 -0.82 2.64 1.00 1455 1724
## y1[8] 5.45 5.45 1.08 1.06 3.73 7.25 1.00 1698 1820
## y1[9] 6.12 6.09 1.02 0.96 4.46 7.82 1.00 1946 1853
## y1[10] 2.69 2.71 1.03 1.03 0.98 4.34 1.00 1851 1676
## y1[11] 8.86 8.82 1.03 1.00 7.17 10.61 1.00 1816 2012
## y1[12] 8.15 8.13 0.98 0.96 6.53 9.77 1.00 2015 1981
## y1[13] 3.86 3.85 1.09 1.03 2.15 5.70 1.00 1831 1859
## y1[14] 9.72 9.75 1.03 0.97 8.04 11.40 1.00 2052 1915
## y1[15] 3.97 3.95 1.04 1.02 2.29 5.66 1.00 1807 1889
## y1[16] 5.93 5.91 1.00 1.01 4.24 7.54 1.00 1941 1801
## y1[17] 0.55 0.53 1.04 0.97 -1.13 2.31 1.00 1378 1802
## y1[18] 7.32 7.30 1.00 1.02 5.68 8.89 1.00 2082 1743
## y1[19] 7.38 7.38 1.03 0.98 5.64 9.05 1.00 2071 1664
## y1[20] 9.94 9.94 1.04 0.99 8.27 11.66 1.00 1898 1952
## y1[21] 9.61 9.58 1.04 1.06 7.94 11.36 1.00 2022 1927
## y1[22] 0.34 0.30 1.11 1.10 -1.41 2.20 1.00 1528 1460
## y1[23] 7.39 7.40 1.02 0.99 5.70 9.00 1.00 1846 1958
## y1[24] 3.12 3.10 1.05 1.03 1.41 4.82 1.00 2090 1968
## y1[25] 4.47 4.47 1.06 1.08 2.72 6.21 1.00 2181 1970
## y1[26] 9.23 9.23 1.03 1.00 7.58 10.90 1.00 1965 1967
## y1[27] 10.19 10.19 1.02 1.00 8.60 11.88 1.00 1957 1735
## y1[28] 1.09 1.06 1.04 1.03 -0.62 2.86 1.00 1559 1734
## y1[29] 10.39 10.40 1.06 1.04 8.64 12.12 1.00 1956 1855
## y1[30] 5.90 5.89 1.03 1.03 4.25 7.60 1.00 1536 1648
## m1[1] 5.15 5.13 0.49 0.50 4.37 6.00 1.00 1114 1261
## m1[2] 4.69 4.69 0.32 0.31 4.17 5.21 1.00 1491 1259
## m1[3] 7.16 7.15 0.44 0.43 6.43 7.89 1.00 972 1186
## m1[4] 6.19 6.17 0.52 0.53 5.37 7.08 1.00 1195 1307
## m1[5] 3.17 3.18 0.33 0.33 2.62 3.70 1.00 760 1026
## m1[6] 5.41 5.41 0.36 0.35 4.81 6.01 1.00 849 1222
## m1[7] 0.93 0.94 0.39 0.38 0.29 1.55 1.00 835 937
## m1[8] 5.44 5.44 0.37 0.35 4.84 6.05 1.00 851 1224
## m1[9] 6.16 6.15 0.26 0.25 5.74 6.57 1.00 1646 1334
## m1[10] 2.64 2.66 0.34 0.33 2.09 3.19 1.00 766 977
## m1[11] 8.84 8.85 0.27 0.26 8.39 9.29 1.00 1709 1267
## m1[12] 8.12 8.12 0.25 0.24 7.68 8.53 1.00 1805 1324
## m1[13] 3.83 3.81 0.48 0.48 3.09 4.65 1.00 1016 1194
## m1[14] 9.72 9.73 0.31 0.29 9.20 10.24 1.00 1581 1274
## m1[15] 3.96 3.96 0.36 0.36 3.38 4.56 1.00 1453 1068
## m1[16] 5.96 5.96 0.26 0.26 5.54 6.39 1.00 1623 1373
## m1[17] 0.59 0.60 0.40 0.39 -0.07 1.24 1.00 855 936
## m1[18] 7.31 7.31 0.24 0.23 6.91 7.72 1.00 1800 1664
## m1[19] 7.37 7.37 0.24 0.23 6.95 7.77 1.00 1805 1608
## m1[20] 9.91 9.92 0.32 0.31 9.37 10.44 1.00 1558 1269
## m1[21] 9.62 9.63 0.31 0.28 9.10 10.12 1.00 1596 1291
## m1[22] 0.29 0.27 0.55 0.56 -0.58 1.21 1.01 905 951
## m1[23] 7.38 7.38 0.24 0.23 6.97 7.79 1.00 1806 1608
## m1[24] 3.12 3.12 0.41 0.41 2.44 3.80 1.00 1402 1028
## m1[25] 4.48 4.48 0.33 0.33 3.94 5.01 1.00 1494 1195
## m1[26] 9.21 9.22 0.29 0.26 8.73 9.68 1.00 1650 1287
## m1[27] 10.17 10.17 0.34 0.32 9.60 10.72 1.00 1533 1297
## m1[28] 1.11 1.13 0.38 0.37 0.49 1.72 1.00 825 928
## m1[29] 10.39 10.39 0.35 0.33 9.80 10.96 1.00 1511 1365
## m1[30] 5.87 5.87 0.38 0.38 5.23 6.50 1.00 880 1361
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ex6
generalized linear regression
poisson regression
objective variable [0,Infinity), integer. variance of error is near to mean
(normal linear regression, correlation is near to 1,-1, variance of error is constant)
# of samples is N,
log li=sum(bj*xji),j=[0,K],i=[1,N]
yi~Po(li)
or
li=sum(bj xji),j=[0,k]
yi~Po(exp li)
when xj -> xj +1, y -> y* exp bj ex6-1.stan
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] l=X*b;
y~poisson_log(l);
}
generated quantities{
array[N] int y1;
vector[N] l1=X*b;
for(i in 1:N){
y1[i]=poisson_log_rng(l1[i]);
}
}n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='poisson')##
## Call: glm(formula = f, family = "poisson", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.347 1.183 0.197
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 52.2
## Residual Deviance: 32 AIC: 96.5X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -71.4455
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -12.1764 0.001982 0.000665165 1 1 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -12.18
## b[1] 0.35
## b[2] 1.18
## b[3] 0.20
## y1[1] 0.00
## y1[2] 1.00
## y1[3] 0.00
## y1[4] 0.00
## y1[5] 5.00
## y1[6] 2.00
## y1[7] 2.00
## y1[8] 1.00
## y1[9] 2.00
## y1[10] 0.00
## y1[11] 3.00
## y1[12] 6.00
## y1[13] 2.00
## y1[14] 0.00
## y1[15] 7.00
## y1[16] 1.00
## y1[17] 4.00
## y1[18] 0.00
## y1[19] 7.00
## y1[20] 6.00
## y1[21] 1.00
## y1[22] 1.00
## y1[23] 2.00
## y1[24] 0.00
## y1[25] 3.00
## y1[26] 2.00
## y1[27] 1.00
## y1[28] 1.00
## y1[29] 2.00
## y1[30] 0.00
## l1[1] -0.07
## l1[2] -0.58
## l1[3] 0.30
## l1[4] -0.06
## l1[5] 1.66
## l1[6] 0.31
## l1[7] 0.78
## l1[8] -0.15
## l1[9] 0.38
## l1[10] -0.15
## l1[11] 1.45
## l1[12] 1.19
## l1[13] -0.24
## l1[14] -0.39
## l1[15] 0.86
## l1[16] -0.16
## l1[17] 0.75
## l1[18] -0.20
## l1[19] 1.12
## l1[20] 1.29
## l1[21] 0.33
## l1[22] 0.31
## l1[23] 1.25
## l1[24] -0.58
## l1[25] 0.86
## l1[26] 0.53
## l1[27] 0.41
## l1[28] 0.29
## l1[29] 0.48
## l1[30] 0.07mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='l1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.65 -13.36 1.20 1.00 -15.91 -12.35 1.01 718 1031
## b[1] 0.32 0.32 0.25 0.25 -0.10 0.72 1.00 967 1092
## b[2] 1.20 1.20 0.27 0.27 0.76 1.65 1.00 1267 1023
## b[3] 0.20 0.20 0.29 0.29 -0.27 0.67 1.00 951 1118
## y1[1] 0.95 1.00 1.01 1.48 0.00 3.00 1.00 2049 2005
## y1[2] 0.55 0.00 0.75 0.00 0.00 2.00 1.00 2024 1959
## y1[3] 1.36 1.00 1.22 1.48 0.00 4.00 1.00 1819 1840
## y1[4] 0.95 1.00 1.00 1.48 0.00 3.00 1.00 1967 1990
## y1[5] 5.24 5.00 2.61 2.97 1.00 10.00 1.00 1873 1910
## y1[6] 1.38 1.00 1.24 1.48 0.00 4.00 1.00 1908 1816
## y1[7] 2.16 2.00 1.60 1.48 0.00 5.00 1.00 1748 1786
## y1[8] 0.90 1.00 0.99 1.48 0.00 3.00 1.00 2033 1940
## y1[9] 1.46 1.00 1.25 1.48 0.00 4.00 1.00 1842 1825
## y1[10] 0.87 1.00 0.97 1.48 0.00 3.00 1.00 2077 2064
## y1[11] 4.34 4.00 2.50 2.97 1.00 9.00 1.00 1608 1335
## y1[12] 3.25 3.00 1.85 1.48 1.00 7.00 1.00 2189 1996
## y1[13] 0.80 1.00 0.94 1.48 0.00 3.00 1.00 1931 1836
## y1[14] 0.68 0.00 0.86 0.00 0.00 2.00 1.00 1800 1768
## y1[15] 2.40 2.00 1.57 1.48 0.00 5.00 1.00 2062 1914
## y1[16] 0.88 1.00 0.97 1.48 0.00 3.00 1.00 1846 1716
## y1[17] 2.08 2.00 1.50 1.48 0.00 5.00 1.00 2095 1971
## y1[18] 0.83 1.00 0.92 1.48 0.00 3.00 1.00 1997 2050
## y1[19] 3.03 3.00 1.84 1.48 0.00 6.00 1.00 2056 2046
## y1[20] 3.70 3.00 1.97 1.48 1.00 7.00 1.00 1904 1896
## y1[21] 1.41 1.00 1.23 1.48 0.00 4.00 1.00 1890 1819
## y1[22] 1.33 1.00 1.18 1.48 0.00 4.00 1.00 2097 1563
## y1[23] 3.49 3.00 1.97 1.48 1.00 7.00 1.00 1881 2004
## y1[24] 0.56 0.00 0.78 0.00 0.00 2.00 1.00 1860 1821
## y1[25] 2.37 2.00 1.60 1.48 0.00 5.00 1.00 1797 1707
## y1[26] 1.65 1.00 1.34 1.48 0.00 4.00 1.00 1883 1938
## y1[27] 1.45 1.00 1.25 1.48 0.00 4.00 1.00 1945 1912
## y1[28] 1.32 1.00 1.19 1.48 0.00 4.00 1.00 1670 1746
## y1[29] 1.65 1.50 1.31 0.74 0.00 4.00 1.00 1917 1910
## y1[30] 1.04 1.00 1.04 1.48 0.00 3.00 1.00 2111 2025
## l1[1] -0.10 -0.10 0.29 0.29 -0.58 0.35 1.00 1007 1220
## l1[2] -0.62 -0.60 0.37 0.38 -1.26 -0.03 1.00 1303 1037
## l1[3] 0.27 0.27 0.21 0.21 -0.10 0.61 1.00 1517 1252
## l1[4] -0.09 -0.08 0.27 0.27 -0.55 0.34 1.00 1380 1142
## l1[5] 1.64 1.65 0.24 0.25 1.23 2.03 1.00 1656 1309
## l1[6] 0.28 0.28 0.25 0.26 -0.15 0.68 1.00 969 1077
## l1[7] 0.76 0.76 0.24 0.25 0.37 1.16 1.00 977 1069
## l1[8] -0.18 -0.18 0.30 0.30 -0.68 0.29 1.00 1018 1196
## l1[9] 0.35 0.35 0.25 0.25 -0.06 0.74 1.00 964 1086
## l1[10] -0.18 -0.17 0.29 0.29 -0.68 0.27 1.00 1358 1033
## l1[11] 1.43 1.43 0.31 0.31 0.92 1.93 1.00 1116 1075
## l1[12] 1.17 1.18 0.18 0.18 0.86 1.44 1.00 1999 1407
## l1[13] -0.27 -0.26 0.30 0.31 -0.80 0.21 1.00 1341 1033
## l1[14] -0.43 -0.42 0.33 0.34 -0.99 0.10 1.00 1051 1189
## l1[15] 0.84 0.85 0.16 0.16 0.57 1.10 1.00 2017 1642
## l1[16] -0.19 -0.19 0.30 0.30 -0.70 0.28 1.00 1020 1166
## l1[17] 0.73 0.73 0.16 0.16 0.45 0.99 1.00 1924 1721
## l1[18] -0.24 -0.23 0.30 0.30 -0.75 0.23 1.00 1348 1044
## l1[19] 1.10 1.11 0.17 0.18 0.81 1.37 1.00 2041 1401
## l1[20] 1.28 1.29 0.19 0.20 0.95 1.57 1.00 1921 1411
## l1[21] 0.30 0.31 0.21 0.20 -0.05 0.64 1.00 1539 1280
## l1[22] 0.28 0.28 0.21 0.21 -0.08 0.62 1.00 1522 1252
## l1[23] 1.23 1.24 0.18 0.19 0.92 1.52 1.00 1954 1432
## l1[24] -0.62 -0.61 0.36 0.36 -1.24 -0.04 1.00 1075 1158
## l1[25] 0.83 0.83 0.25 0.25 0.44 1.24 1.00 990 1101
## l1[26] 0.51 0.51 0.18 0.18 0.20 0.80 1.00 1693 1556
## l1[27] 0.38 0.38 0.24 0.25 -0.03 0.78 1.00 962 1145
## l1[28] 0.26 0.26 0.21 0.21 -0.11 0.60 1.00 1512 1252
## l1[29] 0.46 0.46 0.19 0.19 0.14 0.76 1.00 1651 1536
## l1[30] 0.04 0.05 0.25 0.25 -0.38 0.43 1.00 1417 1174
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 1.5 2 3 4 5
## 0 1 7 0 0 0 0 0
## 1 1 3 0 1 1 0 0
## 2 1 5 0 2 1 0 0
## 3 0 1 1 1 0 0 0
## 4 0 0 0 0 2 0 0
## 6 0 0 0 0 0 1 0
## 7 0 0 0 0 0 0 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 1.5 2 3 4 5
## 0 0 2 0 0 0 0 0
## 1 1 3 0 1 0 0 0
## 2 1 1 0 1 0 0 0
## 3 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0
## 6 0 0 0 0 0 1 0
## 7 0 0 0 0 0 0 0
##
## , , = b
##
##
## 0 1 1.5 2 3 4 5
## 0 1 5 0 0 0 0 0
## 1 0 0 0 0 1 0 0
## 2 0 4 0 1 1 0 0
## 3 0 1 1 1 0 0 0
## 4 0 0 0 0 2 0 0
## 6 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3 5
## 0 4 4 0 0 0
## 1 4 1 1 0 0
## 2 1 6 1 1 0
## 3 1 2 0 0 0
## 4 0 0 0 2 0
## 6 0 0 0 1 0
## 7 0 0 0 0 1cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3 5
## 0 0 2 0 0 0
## 1 4 1 0 0 0
## 2 1 2 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 6 0 0 0 1 0
## 7 0 0 0 0 0
##
## , , = b
##
## map
## 0 1 2 3 5
## 0 4 2 0 0 0
## 1 0 0 1 0 0
## 2 0 4 1 1 0
## 3 1 2 0 0 0
## 4 0 0 0 2 0
## 6 0 0 0 0 0
## 7 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
logistic regression
# of samples is N,
objective variable 0/1 binary
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~Ber(pi), 0/1 binary
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-2.stan
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~bernoulli_logit(z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=bernoulli_rng(inv_logit(z1[i]));
}
}n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='binomial') # it can caluculte when all trials are once##
## Call: glm(formula = f, family = "binomial", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.262 1.222 1.483
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 36.7
## Residual Deviance: 32.6 AIC: 38.6X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -19.4771
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 9 -16.3142 0.000328922 0.000143539 0.954 0.954 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -16.31
## b[1] 0.26
## b[2] 1.22
## b[3] 1.48
## y1[1] 1.00
## y1[2] 0.00
## y1[3] 1.00
## y1[4] 1.00
## y1[5] 1.00
## y1[6] 1.00
## y1[7] 1.00
## y1[8] 0.00
## y1[9] 1.00
## y1[10] 0.00
## y1[11] 0.00
## y1[12] 1.00
## y1[13] 1.00
## y1[14] 1.00
## y1[15] 1.00
## y1[16] 1.00
## y1[17] 1.00
## y1[18] 1.00
## y1[19] 1.00
## y1[20] 1.00
## y1[21] 0.00
## y1[22] 0.00
## y1[23] 0.00
## y1[24] 0.00
## y1[25] 1.00
## y1[26] 1.00
## y1[27] 0.00
## y1[28] 1.00
## y1[29] 1.00
## y1[30] 1.00
## z1[1] 0.47
## z1[2] 0.06
## z1[3] 0.02
## z1[4] 0.83
## z1[5] 1.87
## z1[6] 1.14
## z1[7] 1.46
## z1[8] 1.09
## z1[9] 1.29
## z1[10] 0.00
## z1[11] 0.02
## z1[12] 2.49
## z1[13] 2.36
## z1[14] 1.19
## z1[15] 1.12
## z1[16] 1.86
## z1[17] 1.13
## z1[18] 1.36
## z1[19] -0.62
## z1[20] 2.70
## z1[21] -0.73
## z1[22] 0.57
## z1[23] 0.71
## z1[24] 1.14
## z1[25] 1.01
## z1[26] 0.82
## z1[27] -0.47
## z1[28] 1.43
## z1[29] 0.95
## z1[30] 2.18mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.92 -17.62 1.24 1.05 -20.35 -16.52 1.01 790 1133
## b[1] 0.27 0.27 0.56 0.56 -0.65 1.17 1.00 1536 1568
## b[2] 1.40 1.34 0.92 0.86 -0.04 3.04 1.00 1006 1123
## b[3] 1.76 1.71 1.08 1.08 0.12 3.62 1.01 982 1077
## y1[1] 0.63 1.00 0.48 0.00 0.00 1.00 1.00 1915 NA
## y1[2] 0.53 1.00 0.50 0.00 0.00 1.00 1.00 2065 NA
## y1[3] 0.51 1.00 0.50 0.00 0.00 1.00 1.00 1867 NA
## y1[4] 0.71 1.00 0.45 0.00 0.00 1.00 1.00 2120 NA
## y1[5] 0.87 1.00 0.33 0.00 0.00 1.00 1.00 1874 NA
## y1[6] 0.76 1.00 0.42 0.00 0.00 1.00 1.00 1783 NA
## y1[7] 0.82 1.00 0.38 0.00 0.00 1.00 1.00 2072 NA
## y1[8] 0.75 1.00 0.43 0.00 0.00 1.00 1.00 1966 NA
## y1[9] 0.78 1.00 0.41 0.00 0.00 1.00 1.00 1935 NA
## y1[10] 0.47 0.00 0.50 0.00 0.00 1.00 1.00 1918 NA
## y1[11] 0.50 1.00 0.50 0.00 0.00 1.00 1.00 1995 NA
## y1[12] 0.92 1.00 0.27 0.00 0.00 1.00 1.00 2025 NA
## y1[13] 0.90 1.00 0.30 0.00 0.00 1.00 1.00 1889 NA
## y1[14] 0.77 1.00 0.42 0.00 0.00 1.00 1.00 1900 NA
## y1[15] 0.76 1.00 0.43 0.00 0.00 1.00 1.00 1811 NA
## y1[16] 0.87 1.00 0.34 0.00 0.00 1.00 1.00 1896 NA
## y1[17] 0.76 1.00 0.43 0.00 0.00 1.00 1.00 1923 NA
## y1[18] 0.79 1.00 0.41 0.00 0.00 1.00 1.00 1930 NA
## y1[19] 0.35 0.00 0.48 0.00 0.00 1.00 1.00 1540 NA
## y1[20] 0.93 1.00 0.26 0.00 0.00 1.00 1.00 1914 NA
## y1[21] 0.33 0.00 0.47 0.00 0.00 1.00 1.00 1892 NA
## y1[22] 0.62 1.00 0.49 0.00 0.00 1.00 1.00 1911 NA
## y1[23] 0.67 1.00 0.47 0.00 0.00 1.00 1.00 2071 NA
## y1[24] 0.76 1.00 0.43 0.00 0.00 1.00 1.00 1944 NA
## y1[25] 0.73 1.00 0.44 0.00 0.00 1.00 1.00 2108 NA
## y1[26] 0.71 1.00 0.45 0.00 0.00 1.00 1.00 1882 NA
## y1[27] 0.38 0.00 0.49 0.00 0.00 1.00 1.00 1834 NA
## y1[28] 0.82 1.00 0.38 0.00 0.00 1.00 1.00 2027 NA
## y1[29] 0.73 1.00 0.45 0.00 0.00 1.00 1.00 1847 NA
## y1[30] 0.90 1.00 0.30 0.00 0.00 1.00 1.00 2012 NA
## z1[1] 0.51 0.51 0.55 0.55 -0.38 1.41 1.00 1704 1665
## z1[2] 0.03 0.03 0.61 0.60 -0.98 1.06 1.00 1356 1567
## z1[3] 0.00 -0.01 0.62 0.60 -1.03 1.03 1.00 1329 1542
## z1[4] 0.92 0.90 0.63 0.64 -0.08 1.97 1.00 1727 1401
## z1[5] 2.17 2.08 0.92 0.93 0.80 3.81 1.01 1229 1202
## z1[6] 1.28 1.24 0.77 0.75 0.09 2.59 1.00 1608 1219
## z1[7] 1.70 1.63 0.88 0.88 0.39 3.26 1.00 1380 1188
## z1[8] 1.22 1.19 0.75 0.73 0.08 2.49 1.00 1628 1219
## z1[9] 1.51 1.44 0.89 0.87 0.18 3.08 1.00 1425 1252
## z1[10] -0.03 -0.04 0.63 0.62 -1.09 1.02 1.00 1311 1494
## z1[11] -0.01 -0.01 0.62 0.60 -1.04 1.03 1.00 1328 1518
## z1[12] 2.88 2.79 1.15 1.16 1.15 4.88 1.01 991 1170
## z1[13] 2.73 2.65 1.09 1.11 1.09 4.67 1.01 1014 1216
## z1[14] 1.40 1.34 0.91 0.90 0.04 2.98 1.00 1440 1251
## z1[15] 1.32 1.27 0.92 0.93 -0.06 2.93 1.00 1451 1303
## z1[16] 2.16 2.07 0.92 0.93 0.78 3.78 1.01 1229 1202
## z1[17] 1.26 1.22 0.76 0.74 0.08 2.56 1.00 1616 1219
## z1[18] 1.53 1.48 0.90 0.87 0.17 3.10 1.00 1478 1222
## z1[19] -0.74 -0.73 0.95 0.89 -2.35 0.81 1.00 1027 1238
## z1[20] 3.12 3.01 1.26 1.28 1.22 5.30 1.01 966 1089
## z1[21] -0.86 -0.85 1.02 0.95 -2.63 0.81 1.00 1002 1238
## z1[22] 0.63 0.62 0.56 0.56 -0.29 1.53 1.00 1743 1542
## z1[23] 0.85 0.80 1.06 1.07 -0.81 2.68 1.00 1396 1380
## z1[24] 1.27 1.23 0.77 0.74 0.09 2.58 1.00 1609 1219
## z1[25] 1.12 1.10 0.70 0.70 0.01 2.30 1.00 1662 1305
## z1[26] 0.91 0.89 0.63 0.64 -0.08 1.95 1.00 1730 1401
## z1[27] -0.57 -0.57 0.86 0.82 -2.04 0.86 1.00 1073 1307
## z1[28] 1.68 1.61 0.88 0.87 0.36 3.22 1.00 1387 1188
## z1[29] 1.06 1.03 0.68 0.68 -0.01 2.19 1.00 1684 1286
## z1[30] 2.53 2.44 1.02 1.04 1.00 4.31 1.01 1060 1246
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1
## 0 3 6
## 1 1 20cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1
## 0 3 4
## 1 1 10
##
## , , = b
##
##
## 0 1
## 0 0 2
## 1 0 10par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1
## 0 3 6
## 1 1 20cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1
## 0 3 4
## 1 1 10
##
## , , = b
##
## map
## 0 1
## 0 0 2
## 1 0 10par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
multi logistic regression
# of samples is N,
# of trials of a sample i is mi,
objective variable [0,n], integer
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~B(mi, pi), # of acutual incident
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-3.stan
data{
int N;
int K;
array[N] int m;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~binomial_logit(m,z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
}
}n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)
mdl=cmdstan_model('./ex6-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -96.0634
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -90.0388 0.000161926 5.07823e-05 0.9345 0.9345 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -90.04
## b[1] 0.23
## b[2] 0.60
## b[3] 1.43
## y1[1] 5.00
## y1[2] 8.00
## y1[3] 3.00
## y1[4] 3.00
## y1[5] 7.00
## y1[6] 2.00
## y1[7] 5.00
## y1[8] 2.00
## y1[9] 1.00
## y1[10] 1.00
## y1[11] 2.00
## y1[12] 2.00
## y1[13] 5.00
## y1[14] 1.00
## y1[15] 6.00
## y1[16] 3.00
## y1[17] 3.00
## y1[18] 1.00
## y1[19] 2.00
## y1[20] 8.00
## y1[21] 0.00
## y1[22] 2.00
## y1[23] 4.00
## y1[24] 8.00
## y1[25] 3.00
## y1[26] 2.00
## y1[27] 4.00
## y1[28] 5.00
## y1[29] 4.00
## y1[30] 2.00
## z1[1] 0.61
## z1[2] 2.11
## z1[3] 1.73
## z1[4] 0.68
## z1[5] 2.02
## z1[6] 0.76
## z1[7] 1.91
## z1[8] 0.02
## z1[9] 1.09
## z1[10] 0.12
## z1[11] -0.18
## z1[12] 1.47
## z1[13] 1.57
## z1[14] 0.21
## z1[15] -0.06
## z1[16] 0.18
## z1[17] 1.45
## z1[18] 0.55
## z1[19] 1.60
## z1[20] 2.10
## z1[21] 0.50
## z1[22] -0.21
## z1[23] 1.61
## z1[24] 2.26
## z1[25] 1.45
## z1[26] 0.51
## z1[27] 0.76
## z1[28] 2.01
## z1[29] 0.39
## z1[30] 1.65mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -91.57 -91.27 1.21 1.00 -93.91 -90.23 1.00 901 1189
## b[1] 0.24 0.24 0.23 0.24 -0.13 0.62 1.00 1438 1403
## b[2] 0.59 0.60 0.34 0.35 0.04 1.12 1.00 1491 1325
## b[3] 1.49 1.47 0.41 0.41 0.83 2.16 1.00 1164 910
## y1[1] 5.83 6.00 1.52 1.48 3.00 8.00 1.00 1875 1886
## y1[2] 7.12 7.00 0.93 1.48 5.00 8.00 1.00 2016 NA
## y1[3] 3.42 4.00 0.71 0.00 2.00 4.00 1.00 2021 NA
## y1[4] 2.63 3.00 0.98 1.48 1.00 4.00 1.00 1885 NA
## y1[5] 7.07 7.00 0.95 1.48 5.00 8.00 1.00 1920 NA
## y1[6] 2.72 3.00 0.97 1.48 1.00 4.00 1.00 1969 NA
## y1[7] 5.22 5.00 0.84 1.48 4.00 6.00 1.00 1990 NA
## y1[8] 3.50 4.00 1.41 1.48 1.00 6.00 1.00 1867 1909
## y1[9] 1.50 2.00 0.63 0.00 0.00 2.00 1.00 1898 NA
## y1[10] 1.06 1.00 0.71 1.48 0.00 2.00 1.00 1934 NA
## y1[11] 1.36 1.00 0.89 1.48 0.00 3.00 1.00 1914 NA
## y1[12] 3.29 3.00 0.79 1.48 2.00 4.00 1.00 1811 NA
## y1[13] 6.66 7.00 1.13 1.48 5.00 8.00 1.00 1811 NA
## y1[14] 1.65 2.00 0.88 1.48 0.00 3.00 1.00 1992 NA
## y1[15] 4.37 4.00 1.65 1.48 2.00 7.00 1.00 1851 1938
## y1[16] 3.29 3.00 1.22 1.48 1.00 5.00 1.00 2077 1883
## y1[17] 3.24 3.00 0.81 1.48 2.00 4.00 1.00 1993 NA
## y1[18] 5.04 5.00 1.44 1.48 3.00 7.00 1.00 1994 1945
## y1[19] 1.66 2.00 0.53 0.00 1.00 2.00 1.00 1850 NA
## y1[20] 7.12 7.00 0.92 1.48 5.00 8.00 1.00 1931 NA
## y1[21] 2.50 3.00 1.00 1.48 1.00 4.00 1.00 1629 NA
## y1[22] 3.58 3.50 1.58 2.22 1.00 6.00 1.00 1647 1670
## y1[23] 4.21 4.00 0.85 1.48 3.00 5.00 1.00 1838 NA
## y1[24] 7.24 7.00 0.86 1.48 6.00 8.00 1.00 2019 NA
## y1[25] 2.47 3.00 0.68 0.00 1.00 3.00 1.00 1698 NA
## y1[26] 3.74 4.00 1.20 1.48 2.00 6.00 1.00 1818 NA
## y1[27] 5.42 6.00 1.44 1.48 3.00 8.00 1.00 1901 NA
## y1[28] 4.42 5.00 0.74 0.00 3.00 5.00 1.00 1942 NA
## y1[29] 4.13 4.00 1.35 1.48 2.00 6.00 1.00 1968 1959
## y1[30] 2.53 3.00 0.65 0.00 1.00 3.00 1.00 2003 NA
## z1[1] 0.61 0.61 0.28 0.29 0.15 1.11 1.00 1654 1585
## z1[2] 2.17 2.15 0.38 0.38 1.57 2.82 1.00 1633 1235
## z1[3] 1.79 1.77 0.34 0.33 1.27 2.37 1.00 1453 1346
## z1[4] 0.68 0.68 0.31 0.31 0.17 1.20 1.00 1643 1562
## z1[5] 2.07 2.06 0.36 0.35 1.52 2.68 1.00 1601 1285
## z1[6] 0.76 0.76 0.34 0.35 0.20 1.34 1.00 1635 1563
## z1[7] 1.97 1.95 0.35 0.34 1.43 2.55 1.00 1555 1362
## z1[8] 0.03 0.03 0.28 0.28 -0.43 0.49 1.00 1499 1506
## z1[9] 1.15 1.15 0.52 0.53 0.31 2.02 1.00 1369 1388
## z1[10] 0.13 0.13 0.25 0.26 -0.28 0.54 1.00 1461 1416
## z1[11] -0.17 -0.17 0.35 0.36 -0.76 0.42 1.00 1549 1275
## z1[12] 1.53 1.52 0.38 0.38 0.93 2.16 1.00 1367 1400
## z1[13] 1.63 1.62 0.36 0.36 1.07 2.24 1.00 1392 1346
## z1[14] 0.21 0.21 0.23 0.24 -0.16 0.60 1.00 1442 1393
## z1[15] -0.05 -0.06 0.31 0.31 -0.57 0.46 1.00 1547 1407
## z1[16] 0.18 0.18 0.24 0.25 -0.20 0.57 1.00 1442 1568
## z1[17] 1.51 1.50 0.39 0.39 0.89 2.16 1.00 1366 1444
## z1[18] 0.55 0.55 0.26 0.26 0.12 1.01 1.00 1639 1535
## z1[19] 1.65 1.64 0.36 0.36 1.10 2.25 1.00 1398 1346
## z1[20] 2.15 2.14 0.38 0.38 1.56 2.80 1.00 1631 1251
## z1[21] 0.50 0.49 0.25 0.25 0.09 0.94 1.00 1619 1514
## z1[22] -0.20 -0.20 0.37 0.38 -0.81 0.41 1.00 1546 1275
## z1[23] 1.67 1.65 0.35 0.35 1.12 2.27 1.00 1403 1345
## z1[24] 2.31 2.30 0.43 0.43 1.63 3.05 1.00 1643 1303
## z1[25] 1.51 1.49 0.39 0.39 0.89 2.16 1.00 1365 1422
## z1[26] 0.51 0.51 0.25 0.25 0.10 0.96 1.00 1625 1562
## z1[27] 0.76 0.76 0.34 0.35 0.20 1.34 1.00 1636 1589
## z1[28] 2.06 2.05 0.36 0.35 1.51 2.66 1.00 1598 1235
## z1[29] 0.39 0.39 0.23 0.23 0.02 0.79 1.00 1531 1369
## z1[30] 1.71 1.69 0.35 0.34 1.17 2.30 1.00 1419 1326
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 1 2 3 3.5 4 5 6 7
## 0 0 1 0 0 0 0 0 0
## 1 2 0 0 0 0 0 0 0
## 2 0 1 2 0 2 0 0 0
## 3 0 1 4 0 2 0 0 0
## 4 0 0 2 0 1 1 0 0
## 5 0 0 0 0 1 2 1 1
## 6 0 0 0 1 0 0 0 0
## 7 0 0 0 0 0 0 1 1
## 8 0 0 0 0 0 0 0 3cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 1 2 3 3.5 4 5 6 7
## 0 0 0 0 0 0 0 0 0
## 1 2 0 0 0 0 0 0 0
## 2 0 0 2 0 1 0 0 0
## 3 0 1 1 0 2 0 0 0
## 4 0 0 1 0 0 1 0 0
## 5 0 0 0 0 1 0 1 0
## 6 0 0 0 1 0 0 0 0
## 7 0 0 0 0 0 0 1 0
## 8 0 0 0 0 0 0 0 0
##
## , , = b
##
##
## 1 2 3 3.5 4 5 6 7
## 0 0 1 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0
## 2 0 1 0 0 1 0 0 0
## 3 0 0 3 0 0 0 0 0
## 4 0 0 1 0 1 0 0 0
## 5 0 0 0 0 0 2 0 1
## 6 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 3par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 1 2 3 4 5 6 7 8
## 0 0 1 0 0 0 0 0 0
## 1 2 0 0 0 0 0 0 0
## 2 0 1 2 1 1 0 0 0
## 3 0 1 3 3 0 0 0 0
## 4 0 0 1 2 1 0 0 0
## 5 0 0 0 1 1 2 1 0
## 6 0 0 1 0 0 0 0 0
## 7 0 0 0 0 0 1 0 1
## 8 0 0 0 0 0 0 0 3cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 1 2 3 4 5 6 7 8
## 0 0 0 0 0 0 0 0 0
## 1 2 0 0 0 0 0 0 0
## 2 0 0 2 1 0 0 0 0
## 3 0 1 1 2 0 0 0 0
## 4 0 0 1 0 1 0 0 0
## 5 0 0 0 1 0 1 0 0
## 6 0 0 1 0 0 0 0 0
## 7 0 0 0 0 0 1 0 0
## 8 0 0 0 0 0 0 0 0
##
## , , = b
##
## map
## 1 2 3 4 5 6 7 8
## 0 0 1 0 0 0 0 0 0
## 1 0 0 0 0 0 0 0 0
## 2 0 1 0 0 1 0 0 0
## 3 0 0 2 1 0 0 0 0
## 4 0 0 0 2 0 0 0 0
## 5 0 0 0 0 1 1 1 0
## 6 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 1
## 8 0 0 0 0 0 0 0 3par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
ex7 interaction of variables
n=50
mdl=cmdstan_model('./ex5.stan')
tb=tibble(x=runif(n,-3,3),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,x,y,col=ca),
qplot(data=tb,x,y,col=cb),ncol=2)
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)
fn(f0)## Initial log joint probability = -60.0659
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -20.7022 0.000118055 0.00119508 1 1 19
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -22.83 -22.56 1.44 1.32 -25.52 -21.12 1.00 899 1164
## b[1] 0.14 0.14 0.21 0.20 -0.20 0.49 1.00 607 723
## b[2] 1.22 1.22 0.08 0.08 1.09 1.36 1.00 2429 1901
## b[3] 1.11 1.12 0.29 0.27 0.64 1.58 1.00 493 661
## s 0.97 0.96 0.10 0.10 0.82 1.14 1.00 1622 1486
## y1[1] -0.27 -0.27 0.97 0.95 -1.82 1.38 1.00 1917 1996
## y1[2] -2.41 -2.40 1.04 1.04 -4.13 -0.73 1.00 2205 2024
## y1[3] -1.25 -1.26 0.99 0.99 -2.86 0.39 1.00 2175 2100
## y1[4] -0.71 -0.72 1.00 0.98 -2.39 0.97 1.00 1770 1836
## y1[5] -1.15 -1.16 1.00 0.97 -2.77 0.52 1.00 1776 1757
## y1[6] 2.91 2.90 0.98 0.99 1.31 4.45 1.00 1776 1555
## y1[7] -1.70 -1.67 0.99 1.00 -3.35 -0.11 1.00 1993 1942
## y1[8] 2.90 2.90 0.99 0.94 1.26 4.52 1.00 2028 1933
## y1[9] 1.77 1.76 1.00 0.97 0.14 3.43 1.00 1970 2008
## y1[10] 0.41 0.38 0.98 0.98 -1.19 2.04 1.00 1908 1783
## y1[11] -0.43 -0.43 0.99 0.98 -2.02 1.26 1.00 1947 1882
## y1[12] 2.23 2.24 0.98 0.95 0.62 3.88 1.00 2039 2013
## y1[13] 4.49 4.51 1.06 1.03 2.71 6.18 1.00 1721 1805
## y1[14] 1.93 1.94 1.01 0.97 0.26 3.54 1.00 1924 2010
## y1[15] 0.72 0.70 0.99 0.95 -0.90 2.37 1.00 1775 1801
## y1[16] 0.65 0.64 1.03 1.00 -1.03 2.37 1.00 2033 1831
## y1[17] 3.46 3.48 1.02 0.98 1.78 5.11 1.00 1952 1724
## y1[18] 2.63 2.61 1.00 0.99 1.00 4.30 1.00 1716 1874
## y1[19] -1.02 -1.04 1.00 0.97 -2.66 0.64 1.00 2048 1817
## y1[20] -2.20 -2.15 1.05 1.04 -3.99 -0.50 1.00 2032 1891
## y1[21] -3.29 -3.25 1.02 0.98 -4.98 -1.62 1.00 1811 1852
## y1[22] 2.12 2.14 0.99 0.99 0.48 3.72 1.00 2008 1891
## y1[23] 3.59 3.58 1.01 1.03 1.92 5.25 1.00 1853 1802
## y1[24] 3.51 3.50 0.99 1.00 1.93 5.20 1.00 1944 1952
## y1[25] 2.60 2.59 1.00 0.98 0.96 4.25 1.00 2017 1807
## y1[26] 4.77 4.76 1.02 1.00 3.13 6.41 1.00 1923 1852
## y1[27] -1.44 -1.43 1.01 1.00 -3.08 0.15 1.00 2014 1730
## y1[28] 1.31 1.30 1.01 1.01 -0.35 2.96 1.00 1906 1890
## y1[29] -1.65 -1.68 1.00 1.00 -3.33 0.02 1.00 1891 1892
## y1[30] -1.99 -2.00 0.98 0.97 -3.57 -0.34 1.00 1748 1790
## y1[31] 1.33 1.35 1.03 1.03 -0.34 2.99 1.00 1690 2014
## y1[32] -0.21 -0.21 1.00 0.98 -1.84 1.42 1.00 1754 1593
## y1[33] -0.54 -0.56 0.99 1.01 -2.14 1.02 1.00 2027 1801
## y1[34] -0.54 -0.55 1.01 0.99 -2.23 1.14 1.00 1806 1970
## y1[35] 1.52 1.50 1.02 1.00 -0.14 3.21 1.00 1851 1883
## y1[36] 0.66 0.66 0.98 0.95 -0.91 2.31 1.00 1760 1909
## y1[37] -1.82 -1.82 0.99 0.93 -3.43 -0.18 1.00 2189 1917
## y1[38] 0.87 0.88 1.00 1.01 -0.74 2.51 1.00 1993 1788
## y1[39] 4.27 4.31 1.02 1.01 2.57 5.92 1.00 1866 1953
## y1[40] 3.49 3.51 1.02 1.01 1.71 5.16 1.00 1909 1843
## y1[41] -0.24 -0.22 1.01 0.99 -1.95 1.41 1.00 1914 1970
## y1[42] 2.24 2.24 1.01 1.02 0.64 3.89 1.00 1830 1857
## y1[43] 3.43 3.43 1.01 1.03 1.81 5.08 1.00 1884 1959
## y1[44] 1.16 1.18 1.01 0.94 -0.52 2.82 1.00 1813 1843
## y1[45] -0.75 -0.73 1.01 0.98 -2.43 0.95 1.00 1887 1949
## y1[46] 2.15 2.15 0.98 0.94 0.60 3.74 1.00 1905 1843
## y1[47] 3.60 3.60 1.03 1.06 1.94 5.31 1.00 1949 1939
## y1[48] -1.70 -1.69 0.99 0.96 -3.35 -0.07 1.00 1917 1821
## y1[49] -0.17 -0.18 0.99 0.97 -1.83 1.41 1.00 2093 1969
## y1[50] 2.31 2.30 1.02 1.05 0.64 3.98 1.00 1756 1681
## m1[1] -0.29 -0.29 0.22 0.21 -0.65 0.07 1.00 595 728
## m1[2] -2.40 -2.40 0.29 0.28 -2.90 -1.92 1.00 2403 1543
## m1[3] -1.25 -1.25 0.24 0.23 -1.65 -0.86 1.00 1973 1616
## m1[4] -0.72 -0.72 0.23 0.22 -1.10 -0.34 1.00 600 802
## m1[5] -1.16 -1.16 0.25 0.23 -1.57 -0.76 1.00 618 829
## m1[6] 2.92 2.92 0.23 0.22 2.55 3.30 1.00 1023 1382
## m1[7] -1.69 -1.68 0.26 0.25 -2.12 -1.27 1.00 2161 1618
## m1[8] 2.89 2.88 0.22 0.22 2.52 3.26 1.00 1020 1382
## m1[9] 1.79 1.78 0.22 0.22 1.42 2.15 1.00 858 1027
## m1[10] 0.41 0.41 0.21 0.20 0.08 0.75 1.00 622 703
## m1[11] -0.45 -0.45 0.23 0.22 -0.81 -0.08 1.00 595 776
## m1[12] 2.22 2.21 0.23 0.23 1.84 2.60 1.00 999 1099
## m1[13] 4.45 4.45 0.30 0.29 3.96 4.94 1.00 1208 1554
## m1[14] 1.94 1.93 0.20 0.19 1.61 2.26 1.00 983 1171
## m1[15] 0.68 0.67 0.21 0.20 0.35 1.02 1.00 646 751
## m1[16] 0.65 0.64 0.19 0.18 0.34 0.96 1.00 1185 1340
## m1[17] 3.47 3.47 0.25 0.24 3.06 3.87 1.00 1082 1382
## m1[18] 2.65 2.64 0.25 0.25 2.25 3.04 1.00 1173 1211
## m1[19] -1.03 -1.03 0.23 0.23 -1.42 -0.66 1.00 1879 1638
## m1[20] -2.19 -2.18 0.28 0.27 -2.67 -1.73 1.00 2341 1543
## m1[21] -3.29 -3.28 0.34 0.34 -3.87 -2.72 1.00 793 1095
## m1[22] 2.12 2.11 0.20 0.19 1.79 2.44 1.00 982 1283
## m1[23] 3.60 3.60 0.25 0.25 3.18 4.02 1.00 1096 1445
## m1[24] 3.52 3.51 0.25 0.25 3.10 3.92 1.00 1087 1402
## m1[25] 2.59 2.59 0.24 0.24 2.19 2.98 1.00 1149 1211
## m1[26] 4.80 4.80 0.31 0.31 4.28 5.32 1.00 1240 1597
## m1[27] -1.45 -1.44 0.25 0.24 -1.86 -1.05 1.00 2055 1618
## m1[28] 1.33 1.32 0.21 0.21 0.99 1.67 1.00 746 879
## m1[29] -1.61 -1.60 0.25 0.25 -2.04 -1.19 1.00 2130 1618
## m1[30] -1.98 -1.98 0.28 0.27 -2.45 -1.53 1.00 673 979
## m1[31] 1.32 1.31 0.19 0.18 1.00 1.63 1.00 1034 1179
## m1[32] -0.21 -0.21 0.22 0.21 -0.56 0.15 1.00 596 718
## m1[33] -0.54 -0.54 0.21 0.21 -0.90 -0.19 1.00 1667 1562
## m1[34] -0.49 -0.49 0.23 0.21 -0.86 -0.12 1.00 596 776
## m1[35] 1.51 1.51 0.22 0.21 1.17 1.86 1.00 787 950
## m1[36] 0.68 0.68 0.21 0.20 0.35 1.02 1.00 647 751
## m1[37] -1.85 -1.84 0.27 0.26 -2.30 -1.42 1.00 2224 1645
## m1[38] 0.87 0.87 0.19 0.18 0.57 1.18 1.00 1122 1194
## m1[39] 4.26 4.26 0.29 0.28 3.79 4.73 1.00 1181 1552
## m1[40] 3.50 3.50 0.28 0.27 3.04 3.95 1.00 1555 1645
## m1[41] -0.22 -0.22 0.22 0.21 -0.57 0.14 1.00 595 718
## m1[42] 2.23 2.22 0.23 0.23 1.85 2.61 1.00 1003 1099
## m1[43] 3.42 3.42 0.28 0.27 2.97 3.87 1.00 1518 1587
## m1[44] 1.18 1.18 0.19 0.18 0.87 1.49 1.00 1055 1210
## m1[45] -0.76 -0.76 0.22 0.21 -1.14 -0.40 1.00 1760 1643
## m1[46] 2.13 2.12 0.20 0.19 1.80 2.46 1.00 982 1286
## m1[47] 3.59 3.59 0.25 0.25 3.18 4.01 1.00 1095 1445
## m1[48] -1.71 -1.71 0.26 0.25 -2.15 -1.29 1.00 2169 1618
## m1[49] -0.19 -0.19 0.20 0.20 -0.54 0.14 1.00 1523 1465
## m1[50] 2.32 2.31 0.24 0.23 1.94 2.70 1.00 1036 1098
fn(f1)## Initial log joint probability = -1750.03
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -19.2198 4.63272e-05 0.000150403 1 1 39
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -21.99 -21.62 1.72 1.55 -25.26 -19.91 1.01 797 989
## b[1] -0.07 -0.07 0.25 0.24 -0.50 0.34 1.00 844 972
## b[2] 1.24 1.24 0.08 0.08 1.10 1.37 1.00 2191 1424
## b[3] 1.09 1.09 0.29 0.28 0.64 1.56 1.00 858 871
## b[4] 0.44 0.45 0.27 0.26 -0.02 0.89 1.00 1009 1013
## s 0.95 0.94 0.10 0.10 0.81 1.14 1.00 2295 1549
## y1[1] -0.50 -0.51 0.98 0.96 -2.07 1.08 1.00 1965 1866
## y1[2] -2.65 -2.66 1.03 1.03 -4.32 -0.94 1.00 2034 1685
## y1[3] -1.49 -1.47 1.01 1.00 -3.14 0.14 1.00 1769 1856
## y1[4] -0.51 -0.51 0.97 0.97 -2.14 1.05 1.00 1775 1786
## y1[5] -0.96 -0.97 0.99 0.95 -2.64 0.65 1.00 1707 1761
## y1[6] 2.67 2.69 1.01 1.01 0.99 4.34 1.00 1848 1827
## y1[7] -1.54 -1.55 1.01 1.02 -3.19 0.13 1.00 1932 1936
## y1[8] 3.12 3.11 0.98 0.97 1.50 4.76 1.00 1905 1971
## y1[9] 2.04 2.06 1.02 1.01 0.38 3.71 1.00 1879 1925
## y1[10] 0.62 0.63 0.97 0.98 -0.94 2.20 1.00 1930 1989
## y1[11] -0.67 -0.69 1.00 0.99 -2.33 1.01 1.00 1878 1730
## y1[12] 2.03 2.04 1.02 1.00 0.31 3.66 1.00 1996 2093
## y1[13] 4.66 4.66 1.03 1.01 2.97 6.35 1.00 1959 1973
## y1[14] 1.66 1.65 0.98 0.98 0.05 3.24 1.00 1940 1972
## y1[15] 0.92 0.94 0.99 1.00 -0.72 2.48 1.00 1900 1890
## y1[16] 0.85 0.86 0.97 0.99 -0.77 2.41 1.00 2094 1743
## y1[17] 3.69 3.70 1.02 0.99 2.06 5.39 1.00 1907 1828
## y1[18] 2.91 2.88 1.02 1.01 1.27 4.59 1.00 2080 2059
## y1[19] -1.27 -1.26 1.01 1.02 -2.94 0.35 1.00 1952 1945
## y1[20] -2.00 -1.97 1.00 0.97 -3.66 -0.41 1.00 1898 1813
## y1[21] -3.51 -3.53 1.02 0.99 -5.20 -1.87 1.00 1746 1783
## y1[22] 1.84 1.84 0.99 0.99 0.19 3.47 1.00 1919 1858
## y1[23] 3.38 3.37 1.00 0.96 1.76 5.10 1.00 1966 1852
## y1[24] 3.29 3.31 1.01 0.98 1.63 4.93 1.00 1968 1801
## y1[25] 2.41 2.42 1.00 0.97 0.72 4.09 1.00 1737 1931
## y1[26] 5.09 5.12 1.03 1.00 3.39 6.76 1.00 1989 1984
## y1[27] -1.26 -1.28 1.01 1.00 -2.86 0.43 1.00 1955 1856
## y1[28] 1.11 1.11 0.99 0.97 -0.57 2.72 1.00 1856 1877
## y1[29] -1.49 -1.48 1.00 1.01 -3.14 0.15 1.00 1907 1832
## y1[30] -1.75 -1.79 1.00 0.99 -3.34 -0.06 1.00 1806 1483
## y1[31] 1.09 1.10 0.96 0.94 -0.47 2.66 1.00 1883 1907
## y1[32] -0.41 -0.42 0.98 0.95 -1.99 1.22 1.00 1683 1816
## y1[33] -0.80 -0.77 1.00 0.97 -2.48 0.84 1.00 1987 1740
## y1[34] -0.71 -0.73 0.99 0.98 -2.29 0.95 1.00 1876 1927
## y1[35] 1.29 1.28 0.96 0.94 -0.23 2.90 1.00 1578 1909
## y1[36] 0.92 0.91 1.00 0.97 -0.72 2.53 1.00 1910 1970
## y1[37] -1.68 -1.68 0.98 0.99 -3.35 -0.03 1.00 2019 2010
## y1[38] 0.61 0.60 1.00 0.97 -1.03 2.22 1.00 1917 1865
## y1[39] 4.05 4.06 1.00 0.98 2.41 5.70 1.00 2006 2084
## y1[40] 3.83 3.84 1.02 0.99 2.12 5.48 1.00 1732 1854
## y1[41] -0.42 -0.45 0.99 0.94 -2.00 1.19 1.00 2013 1806
## y1[42] 2.49 2.50 0.99 0.96 0.83 4.11 1.00 1686 1934
## y1[43] 3.26 3.27 0.99 0.99 1.65 4.87 1.00 1807 1850
## y1[44] 1.39 1.39 0.98 0.97 -0.20 3.04 1.00 2083 1965
## y1[45] -0.56 -0.56 0.97 0.96 -2.22 1.01 1.00 1930 1972
## y1[46] 2.35 2.34 1.01 1.00 0.61 4.01 1.00 1899 1902
## y1[47] 3.38 3.37 0.98 0.94 1.73 4.94 1.00 1682 1810
## y1[48] -1.53 -1.53 0.98 0.96 -3.14 0.08 1.00 1915 1815
## y1[49] -0.01 -0.04 0.99 0.98 -1.62 1.64 1.00 2119 1844
## y1[50] 2.12 2.13 0.98 1.00 0.56 3.70 1.00 1900 1774
## m1[1] -0.51 -0.51 0.26 0.25 -0.95 -0.10 1.00 845 953
## m1[2] -2.68 -2.68 0.34 0.34 -3.25 -2.12 1.00 1846 1847
## m1[3] -1.51 -1.51 0.29 0.28 -1.99 -1.03 1.00 1613 1473
## m1[4] -0.50 -0.49 0.26 0.25 -0.93 -0.07 1.00 957 1180
## m1[5] -0.94 -0.95 0.27 0.26 -1.38 -0.50 1.00 966 1187
## m1[6] 2.70 2.70 0.26 0.25 2.28 3.12 1.00 1119 1313
## m1[7] -1.51 -1.51 0.28 0.27 -1.98 -1.07 1.00 2459 1807
## m1[8] 3.11 3.11 0.27 0.26 2.66 3.54 1.00 1568 1294
## m1[9] 2.03 2.03 0.27 0.26 1.60 2.46 1.00 1246 1590
## m1[10] 0.65 0.64 0.25 0.24 0.22 1.06 1.00 1024 1254
## m1[11] -0.67 -0.66 0.26 0.26 -1.12 -0.25 1.00 850 950
## m1[12] 2.02 2.02 0.26 0.26 1.59 2.47 1.00 1028 1271
## m1[13] 4.69 4.70 0.34 0.33 4.13 5.24 1.00 1637 1461
## m1[14] 1.70 1.71 0.24 0.24 1.31 2.08 1.00 1107 1177
## m1[15] 0.91 0.91 0.25 0.24 0.49 1.33 1.00 1054 1346
## m1[16] 0.85 0.85 0.22 0.21 0.48 1.20 1.00 1814 1416
## m1[17] 3.70 3.70 0.29 0.28 3.20 4.17 1.00 1584 1327
## m1[18] 2.90 2.90 0.29 0.28 2.43 3.39 1.00 1433 1616
## m1[19] -1.30 -1.29 0.28 0.28 -1.75 -0.83 1.00 1569 1526
## m1[20] -2.02 -2.01 0.30 0.30 -2.52 -1.55 1.00 2568 1874
## m1[21] -3.54 -3.53 0.37 0.37 -4.16 -2.93 1.00 1077 1382
## m1[22] 1.89 1.89 0.24 0.24 1.49 2.27 1.00 1103 1176
## m1[23] 3.38 3.39 0.29 0.28 2.92 3.85 1.00 1151 1482
## m1[24] 3.30 3.31 0.28 0.27 2.84 3.76 1.00 1146 1482
## m1[25] 2.40 2.40 0.27 0.27 1.94 2.85 1.00 1094 1277
## m1[26] 5.04 5.05 0.35 0.35 4.45 5.62 1.00 1659 1525
## m1[27] -1.27 -1.26 0.27 0.26 -1.71 -0.84 1.00 2406 1794
## m1[28] 1.13 1.13 0.25 0.24 0.70 1.54 1.00 910 1100
## m1[29] -1.43 -1.43 0.27 0.27 -1.89 -1.00 1.00 2441 1841
## m1[30] -1.77 -1.77 0.30 0.29 -2.26 -1.30 1.00 999 1144
## m1[31] 1.08 1.08 0.24 0.23 0.69 1.46 1.00 1150 1296
## m1[32] -0.43 -0.42 0.26 0.25 -0.86 -0.02 1.00 844 961
## m1[33] -0.80 -0.79 0.26 0.26 -1.23 -0.37 1.00 1474 1491
## m1[34] -0.71 -0.71 0.26 0.26 -1.17 -0.29 1.00 852 970
## m1[35] 1.31 1.31 0.25 0.25 0.88 1.73 1.00 930 1139
## m1[36] 0.92 0.92 0.25 0.24 0.50 1.33 1.00 1054 1346
## m1[37] -1.68 -1.67 0.28 0.28 -2.15 -1.23 1.00 2495 1731
## m1[38] 0.63 0.63 0.24 0.23 0.24 1.01 1.00 1196 1277
## m1[39] 4.06 4.06 0.31 0.31 3.54 4.56 1.00 1193 1533
## m1[40] 3.76 3.76 0.33 0.32 3.22 4.31 1.00 1631 1711
## m1[41] -0.44 -0.43 0.26 0.25 -0.87 -0.03 1.00 844 953
## m1[42] 2.48 2.48 0.28 0.27 2.03 2.94 1.00 1340 1623
## m1[43] 3.24 3.23 0.30 0.30 2.73 3.75 1.00 1265 1324
## m1[44] 1.39 1.39 0.23 0.21 1.02 1.75 1.00 1694 1386
## m1[45] -0.58 -0.57 0.24 0.24 -0.98 -0.18 1.00 2223 1639
## m1[46] 2.34 2.34 0.24 0.24 1.94 2.74 1.00 1587 1350
## m1[47] 3.38 3.38 0.29 0.28 2.91 3.84 1.00 1151 1482
## m1[48] -1.54 -1.53 0.28 0.27 -2.00 -1.10 1.00 2464 1790
## m1[49] 0.00 0.00 0.23 0.23 -0.39 0.37 1.00 2050 1618
## m1[50] 2.12 2.12 0.27 0.26 1.68 2.57 1.00 1047 1271
fn(f2)## Initial log joint probability = -287.166
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 36 -18.3444 5.62153e-05 0.00205695 0.9495 0.9495 42
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -21.67 -21.38 1.81 1.69 -25.02 -19.37 1.01 804 1047
## b[1] -0.24 -0.25 0.27 0.26 -0.69 0.20 1.01 642 1004
## b[2] 1.23 1.23 0.08 0.08 1.09 1.37 1.00 2428 1425
## b[3] 1.43 1.45 0.37 0.36 0.83 2.01 1.01 527 951
## b[4] 0.81 0.81 0.40 0.39 0.18 1.48 1.01 558 904
## b[5] -0.68 -0.68 0.53 0.54 -1.57 0.19 1.02 428 753
## s 0.95 0.94 0.10 0.10 0.80 1.12 1.00 1989 1462
## y1[1] -0.69 -0.69 1.02 0.97 -2.41 0.94 1.00 1848 1954
## y1[2] -2.47 -2.48 1.05 1.03 -4.19 -0.78 1.00 1730 2029
## y1[3] -1.29 -1.28 0.98 0.96 -2.89 0.34 1.00 1875 1767
## y1[4] -0.31 -0.30 0.97 0.98 -1.94 1.25 1.00 2050 1929
## y1[5] -0.74 -0.72 1.00 0.98 -2.41 0.82 1.00 1938 2012
## y1[6] 2.84 2.86 1.02 1.02 1.20 4.52 1.00 1910 1678
## y1[7] -1.60 -1.59 0.99 0.97 -3.26 0.06 1.00 2020 1888
## y1[8] 2.97 3.00 0.99 0.97 1.32 4.54 1.00 1782 2003
## y1[9] 2.22 2.20 0.99 0.98 0.56 3.86 1.00 1951 1820
## y1[10] 0.81 0.78 0.98 0.96 -0.78 2.45 1.00 2198 2008
## y1[11] -0.83 -0.82 1.00 0.97 -2.49 0.78 1.00 2129 2148
## y1[12] 1.82 1.81 1.00 0.97 0.22 3.47 1.00 1823 1810
## y1[13] 4.52 4.54 1.02 1.05 2.82 6.21 1.00 1728 1804
## y1[14] 1.86 1.88 0.97 0.93 0.31 3.44 1.00 1926 1995
## y1[15] 1.09 1.05 0.99 0.97 -0.51 2.73 1.00 1844 1945
## y1[16] 0.74 0.74 0.97 0.95 -0.84 2.39 1.00 2092 1905
## y1[17] 3.55 3.57 1.00 0.97 1.90 5.23 1.00 1972 1918
## y1[18] 3.12 3.12 1.00 1.00 1.48 4.77 1.00 1886 1815
## y1[19] -1.13 -1.11 1.03 0.99 -2.90 0.55 1.00 1994 1774
## y1[20] -2.15 -2.14 1.01 0.98 -3.83 -0.45 1.00 2118 1680
## y1[21] -3.71 -3.74 1.04 1.04 -5.36 -1.98 1.00 1601 1902
## y1[22] 2.04 2.03 0.96 0.94 0.51 3.66 1.00 1905 1810
## y1[23] 3.54 3.55 1.00 1.01 1.90 5.14 1.00 1943 1857
## y1[24] 3.46 3.45 1.01 1.01 1.80 5.16 1.00 2011 1834
## y1[25] 2.23 2.20 0.98 0.96 0.68 3.88 1.00 1925 2040
## y1[26] 4.83 4.82 1.01 1.04 3.16 6.53 1.00 1968 1857
## y1[27] -1.39 -1.38 0.99 0.97 -3.05 0.24 1.00 2144 1878
## y1[28] 0.94 0.95 1.01 0.99 -0.75 2.55 1.00 1945 1968
## y1[29] -1.50 -1.50 0.98 0.98 -3.09 0.12 1.00 1978 1934
## y1[30] -1.56 -1.56 1.02 1.00 -3.21 0.13 1.00 1847 1874
## y1[31] 1.22 1.25 0.99 0.95 -0.46 2.77 1.00 1839 1911
## y1[32] -0.58 -0.56 1.00 1.01 -2.21 1.04 1.00 1895 1811
## y1[33] -0.59 -0.61 1.00 0.97 -2.23 1.03 1.00 1745 1883
## y1[34] -0.90 -0.90 0.97 0.94 -2.49 0.71 1.00 1987 1892
## y1[35] 1.11 1.13 0.96 0.96 -0.50 2.64 1.00 1654 1929
## y1[36] 1.13 1.15 0.97 0.96 -0.46 2.68 1.00 1918 2081
## y1[37] -1.79 -1.76 0.99 0.97 -3.45 -0.19 1.00 1917 1932
## y1[38] 0.81 0.79 0.98 0.96 -0.82 2.41 1.00 1770 1767
## y1[39] 4.20 4.19 0.99 0.98 2.58 5.81 1.00 2060 1675
## y1[40] 3.95 3.94 1.03 1.02 2.27 5.69 1.00 1966 1629
## y1[41] -0.62 -0.63 0.99 0.97 -2.28 1.05 1.00 1886 1853
## y1[42] 2.67 2.68 1.01 1.03 0.98 4.31 1.00 2018 1931
## y1[43] 3.02 3.01 0.99 0.96 1.36 4.62 1.00 1945 2014
## y1[44] 1.24 1.25 0.98 0.97 -0.39 2.84 1.00 1955 1906
## y1[45] -0.68 -0.67 0.94 0.93 -2.26 0.82 1.00 2112 1880
## y1[46] 2.20 2.19 1.00 0.96 0.51 3.83 1.00 1894 1866
## y1[47] 3.53 3.52 1.02 1.00 1.86 5.20 1.00 2037 1932
## y1[48] -1.66 -1.66 0.98 0.98 -3.24 0.00 1.00 1837 1730
## y1[49] -0.11 -0.09 0.99 0.97 -1.81 1.46 1.00 2107 2015
## y1[50] 1.90 1.91 0.99 1.00 0.31 3.55 1.00 1966 1890
## m1[1] -0.68 -0.68 0.27 0.28 -1.14 -0.22 1.01 649 1003
## m1[2] -2.47 -2.47 0.38 0.37 -3.07 -1.84 1.00 1440 1542
## m1[3] -1.32 -1.32 0.33 0.33 -1.84 -0.78 1.00 1315 1329
## m1[4] -0.29 -0.29 0.31 0.31 -0.81 0.23 1.00 1176 1337
## m1[5] -0.73 -0.73 0.32 0.32 -1.26 -0.20 1.00 1208 1396
## m1[6] 2.86 2.86 0.28 0.28 2.40 3.31 1.01 1209 1398
## m1[7] -1.62 -1.63 0.30 0.28 -2.11 -1.11 1.00 1896 1398
## m1[8] 2.96 2.95 0.30 0.31 2.49 3.44 1.00 1580 1477
## m1[9] 2.22 2.22 0.30 0.30 1.71 2.73 1.01 1172 1295
## m1[10] 0.85 0.84 0.30 0.29 0.35 1.35 1.01 1149 1177
## m1[11] -0.83 -0.83 0.28 0.28 -1.29 -0.37 1.00 655 1002
## m1[12] 1.84 1.84 0.29 0.28 1.36 2.32 1.00 745 1277
## m1[13] 4.53 4.52 0.37 0.39 3.94 5.13 1.00 1683 1539
## m1[14] 1.87 1.88 0.26 0.26 1.43 2.30 1.01 1140 1403
## m1[15] 1.11 1.11 0.29 0.28 0.62 1.61 1.01 1144 1271
## m1[16] 0.71 0.71 0.26 0.25 0.31 1.16 1.00 1568 1497
## m1[17] 3.54 3.53 0.33 0.34 3.02 4.07 1.00 1617 1434
## m1[18] 3.08 3.07 0.32 0.32 2.54 3.62 1.01 1241 1308
## m1[19] -1.10 -1.10 0.32 0.32 -1.61 -0.58 1.00 1290 1316
## m1[20] -2.12 -2.13 0.32 0.30 -2.64 -1.58 1.00 1995 1315
## m1[21] -3.67 -3.67 0.38 0.38 -4.30 -3.04 1.00 897 1466
## m1[22] 2.05 2.05 0.27 0.27 1.62 2.48 1.01 1149 1417
## m1[23] 3.53 3.54 0.30 0.31 3.05 4.02 1.01 1302 1380
## m1[24] 3.45 3.45 0.30 0.31 2.96 3.93 1.01 1291 1346
## m1[25] 2.21 2.21 0.30 0.29 1.71 2.71 1.00 786 1220
## m1[26] 4.87 4.86 0.39 0.40 4.25 5.51 1.00 1703 1491
## m1[27] -1.38 -1.39 0.29 0.27 -1.85 -0.89 1.00 1852 1384
## m1[28] 0.94 0.94 0.27 0.26 0.50 1.41 1.00 673 1058
## m1[29] -1.54 -1.55 0.29 0.27 -2.02 -1.04 1.00 1886 1406
## m1[30] -1.55 -1.55 0.35 0.35 -2.13 -0.98 1.00 1283 1512
## m1[31] 1.25 1.25 0.26 0.27 0.81 1.68 1.01 1127 1352
## m1[32] -0.59 -0.60 0.27 0.27 -1.05 -0.14 1.01 647 972
## m1[33] -0.61 -0.61 0.30 0.29 -1.09 -0.11 1.01 1236 1304
## m1[34] -0.87 -0.88 0.28 0.28 -1.34 -0.42 1.00 655 1003
## m1[35] 1.13 1.13 0.27 0.27 0.69 1.60 1.00 685 1167
## m1[36] 1.11 1.11 0.29 0.28 0.62 1.61 1.01 1144 1271
## m1[37] -1.78 -1.79 0.30 0.28 -2.28 -1.26 1.00 1913 1419
## m1[38] 0.80 0.81 0.27 0.26 0.36 1.24 1.01 1135 1439
## m1[39] 4.20 4.20 0.32 0.33 3.67 4.72 1.00 1417 1537
## m1[40] 3.93 3.93 0.35 0.35 3.35 4.51 1.01 1329 1340
## m1[41] -0.60 -0.61 0.27 0.27 -1.06 -0.15 1.01 647 972
## m1[42] 2.66 2.66 0.31 0.31 2.14 3.17 1.01 1207 1265
## m1[43] 3.04 3.03 0.33 0.32 2.50 3.59 1.00 897 1296
## m1[44] 1.25 1.25 0.26 0.25 0.83 1.69 1.00 1519 1528
## m1[45] -0.69 -0.70 0.27 0.25 -1.13 -0.23 1.00 1762 1497
## m1[46] 2.20 2.20 0.28 0.29 1.75 2.64 1.00 1537 1595
## m1[47] 3.53 3.53 0.30 0.31 3.04 4.02 1.01 1301 1380
## m1[48] -1.65 -1.66 0.30 0.28 -2.14 -1.14 1.00 1900 1420
## m1[49] -0.13 -0.13 0.26 0.25 -0.54 0.32 1.00 1678 1547
## m1[50] 1.93 1.93 0.29 0.28 1.45 2.42 1.00 754 1254
tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,xa,y,col=ca),
qplot(data=tb,xa,y,col=cb),
qplot(data=tb,xb,y,col=ca),
qplot(data=tb,xb,y,col=cb))
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)
fn(f0)## Initial log joint probability = -1547.39
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 -45.4497 0.00023767 0.000932955 1 1 30
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -48.25 -47.89 1.85 1.67 -51.79 -45.94 1.01 642 1094
## b[1] 0.34 0.35 0.41 0.40 -0.33 1.01 1.01 913 999
## b[2] 1.09 1.09 0.20 0.19 0.78 1.41 1.00 2680 1372
## b[3] 1.14 1.14 0.24 0.24 0.75 1.53 1.00 2296 1334
## b[4] 0.77 0.77 0.49 0.47 -0.06 1.53 1.01 877 1010
## b[5] 1.65 1.64 0.48 0.45 0.88 2.50 1.01 800 959
## s 1.63 1.61 0.18 0.17 1.36 1.95 1.00 1928 1676
## y1[1] 0.09 0.11 1.70 1.71 -2.64 2.88 1.00 1754 1763
## y1[2] 1.14 1.12 1.71 1.73 -1.57 3.86 1.00 1855 1683
## y1[3] 0.14 0.16 1.69 1.66 -2.66 2.85 1.00 1917 2104
## y1[4] 5.95 5.93 1.69 1.63 3.22 8.72 1.00 2023 1573
## y1[5] -1.87 -1.88 1.66 1.64 -4.59 0.83 1.00 1836 1890
## y1[6] 3.29 3.31 1.71 1.67 0.46 6.05 1.00 1886 1681
## y1[7] 4.34 4.37 1.72 1.69 1.51 7.25 1.00 2013 1857
## y1[8] 2.84 2.85 1.70 1.66 0.07 5.58 1.00 2126 2102
## y1[9] -0.70 -0.72 1.73 1.72 -3.58 2.17 1.00 1944 1742
## y1[10] 0.51 0.56 1.70 1.71 -2.33 3.19 1.00 1851 1916
## y1[11] 2.82 2.86 1.68 1.67 0.05 5.52 1.00 1889 1582
## y1[12] 4.93 4.93 1.71 1.68 2.16 7.71 1.00 2157 1853
## y1[13] 0.23 0.23 1.72 1.72 -2.56 3.11 1.00 1997 1974
## y1[14] 0.81 0.84 1.74 1.71 -2.02 3.65 1.00 2004 2053
## y1[15] -2.39 -2.42 1.78 1.74 -5.25 0.53 1.00 1957 1873
## y1[16] -1.16 -1.13 1.72 1.74 -4.01 1.66 1.00 1954 1989
## y1[17] 3.55 3.56 1.74 1.71 0.57 6.34 1.00 1989 1784
## y1[18] 1.37 1.31 1.67 1.67 -1.30 4.23 1.00 2001 1937
## y1[19] 2.66 2.62 1.76 1.70 -0.20 5.60 1.00 2107 1856
## y1[20] 3.61 3.60 1.71 1.72 0.84 6.48 1.00 2033 1930
## y1[21] -0.40 -0.35 1.77 1.74 -3.30 2.47 1.00 1889 2004
## y1[22] 3.12 3.14 1.71 1.69 0.28 5.94 1.00 1675 1875
## y1[23] 3.39 3.40 1.67 1.64 0.55 6.19 1.00 1932 1890
## y1[24] 2.90 2.86 1.77 1.74 0.08 5.78 1.00 1648 1818
## y1[25] 1.15 1.13 1.78 1.77 -1.74 4.07 1.00 1854 1686
## y1[26] -1.56 -1.55 1.73 1.65 -4.43 1.31 1.00 1969 1895
## y1[27] 3.97 3.98 1.72 1.70 1.20 6.78 1.00 1897 1906
## y1[28] 1.41 1.40 1.73 1.69 -1.41 4.32 1.00 1913 1855
## y1[29] 1.22 1.21 1.71 1.67 -1.57 4.13 1.00 1811 1962
## y1[30] -2.31 -2.36 1.76 1.74 -5.10 0.54 1.00 1725 1523
## y1[31] 1.96 2.00 1.72 1.70 -0.85 4.71 1.00 1818 1843
## y1[32] 0.20 0.19 1.73 1.71 -2.62 3.02 1.00 1918 1934
## y1[33] 0.41 0.42 1.73 1.78 -2.37 3.24 1.00 2052 1973
## y1[34] 2.16 2.15 1.73 1.70 -0.65 5.09 1.00 1839 1921
## y1[35] -1.24 -1.24 1.67 1.67 -4.02 1.52 1.00 2028 1844
## y1[36] 1.48 1.47 1.76 1.70 -1.39 4.36 1.00 1908 2004
## y1[37] 1.10 1.07 1.67 1.69 -1.61 3.83 1.00 1671 1792
## y1[38] -1.20 -1.20 1.72 1.66 -4.15 1.53 1.00 1623 1807
## y1[39] 2.09 2.07 1.67 1.65 -0.72 4.81 1.00 1983 1884
## y1[40] 3.06 3.05 1.75 1.72 0.21 5.89 1.00 2188 1961
## y1[41] -1.94 -1.92 1.74 1.68 -4.88 0.80 1.00 2201 1991
## y1[42] 3.77 3.82 1.70 1.64 0.91 6.53 1.00 1904 1678
## y1[43] 1.97 2.01 1.72 1.67 -0.85 4.79 1.00 1951 1970
## y1[44] -0.94 -0.91 1.67 1.69 -3.69 1.76 1.00 2202 1917
## y1[45] -0.05 -0.04 1.65 1.64 -2.75 2.66 1.00 1912 1856
## y1[46] -0.04 -0.04 1.75 1.69 -2.88 2.82 1.00 1916 1704
## y1[47] 3.54 3.54 1.67 1.66 0.75 6.20 1.00 2159 1933
## y1[48] -2.94 -2.95 1.74 1.76 -5.82 -0.06 1.00 1613 1874
## y1[49] -1.24 -1.23 1.71 1.67 -4.00 1.56 1.00 2070 1917
## y1[50] 0.82 0.82 1.67 1.57 -1.94 3.53 1.00 1790 1969
## m1[1] 0.09 0.09 0.49 0.50 -0.75 0.86 1.00 1370 1308
## m1[2] 1.17 1.20 0.38 0.39 0.52 1.76 1.00 1219 1314
## m1[3] 0.13 0.15 0.38 0.38 -0.50 0.74 1.00 1264 1411
## m1[4] 5.96 5.95 0.59 0.59 5.01 6.93 1.01 2195 1494
## m1[5] -1.82 -1.83 0.49 0.47 -2.63 -1.03 1.00 2057 1437
## m1[6] 3.33 3.34 0.40 0.40 2.68 4.00 1.00 1651 1463
## m1[7] 4.37 4.38 0.67 0.66 3.26 5.46 1.00 1661 1608
## m1[8] 2.84 2.84 0.54 0.53 1.97 3.71 1.00 1342 1430
## m1[9] -0.67 -0.68 0.51 0.51 -1.49 0.17 1.00 1290 1372
## m1[10] 0.51 0.51 0.49 0.49 -0.29 1.29 1.00 1508 1634
## m1[11] 2.79 2.81 0.47 0.48 2.01 3.53 1.00 1309 1537
## m1[12] 4.90 4.90 0.51 0.50 4.06 5.74 1.01 1795 1469
## m1[13] 0.24 0.24 0.47 0.47 -0.53 0.99 1.00 2029 1595
## m1[14] 0.90 0.89 0.67 0.65 -0.17 2.02 1.00 1899 1401
## m1[15] -2.41 -2.40 0.53 0.51 -3.31 -1.57 1.00 2032 1458
## m1[16] -1.17 -1.15 0.45 0.45 -1.90 -0.45 1.00 1216 1281
## m1[17] 3.57 3.58 0.53 0.55 2.68 4.40 1.00 1525 1645
## m1[18] 1.34 1.35 0.41 0.42 0.67 1.99 1.00 1729 1336
## m1[19] 2.65 2.64 0.64 0.63 1.61 3.71 1.00 1691 1382
## m1[20] 3.56 3.57 0.50 0.49 2.75 4.39 1.00 1418 1459
## m1[21] -0.33 -0.34 0.64 0.63 -1.39 0.71 1.00 1899 1524
## m1[22] 3.06 3.05 0.52 0.50 2.23 3.92 1.01 1258 1109
## m1[23] 3.39 3.40 0.41 0.41 2.72 4.06 1.00 1511 1348
## m1[24] 2.91 2.91 0.58 0.57 1.96 3.85 1.00 1219 1108
## m1[25] 1.17 1.15 0.64 0.63 0.11 2.21 1.00 1423 1027
## m1[26] -1.51 -1.50 0.46 0.44 -2.27 -0.77 1.00 1681 1551
## m1[27] 3.96 3.97 0.53 0.54 3.13 4.86 1.00 2621 1646
## m1[28] 1.43 1.46 0.40 0.40 0.77 2.06 1.00 1333 1361
## m1[29] 1.25 1.25 0.61 0.60 0.22 2.22 1.01 1115 1169
## m1[30] -2.31 -2.30 0.54 0.54 -3.19 -1.45 1.00 1186 1402
## m1[31] 1.97 1.96 0.61 0.62 1.00 2.98 1.01 1155 1289
## m1[32] 0.19 0.19 0.52 0.51 -0.66 1.03 1.00 1533 1662
## m1[33] 0.47 0.45 0.53 0.52 -0.39 1.35 1.01 1219 1290
## m1[34] 2.13 2.15 0.53 0.54 1.22 2.99 1.01 1121 1421
## m1[35] -1.28 -1.27 0.64 0.62 -2.36 -0.19 1.00 1459 1416
## m1[36] 1.43 1.42 0.58 0.59 0.49 2.36 1.01 1456 1343
## m1[37] 1.12 1.12 0.43 0.44 0.42 1.80 1.00 1483 1402
## m1[38] -1.26 -1.26 0.49 0.48 -2.07 -0.48 1.00 2176 1363
## m1[39] 2.03 2.06 0.42 0.42 1.33 2.68 1.00 1187 1365
## m1[40] 3.08 3.08 0.51 0.52 2.23 3.89 1.00 1366 1302
## m1[41] -1.91 -1.91 0.57 0.57 -2.87 -0.95 1.00 1572 1539
## m1[42] 3.75 3.74 0.51 0.52 2.95 4.61 1.00 2544 1519
## m1[43] 1.95 1.94 0.58 0.59 1.00 2.94 1.00 1635 1284
## m1[44] -0.97 -0.96 0.43 0.41 -1.68 -0.28 1.00 1709 1433
## m1[45] -0.07 -0.07 0.41 0.40 -0.75 0.60 1.01 901 975
## m1[46] -0.06 -0.07 0.59 0.58 -1.01 0.93 1.00 1414 1381
## m1[47] 3.56 3.57 0.54 0.56 2.71 4.46 1.00 2650 1653
## m1[48] -2.95 -2.95 0.58 0.58 -3.87 -2.05 1.00 1339 1536
## m1[49] -1.21 -1.20 0.44 0.42 -1.93 -0.50 1.00 1725 1543
## m1[50] 0.84 0.86 0.39 0.39 0.19 1.45 1.01 1089 1400
fn(f1)## Initial log joint probability = -100.856
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 36 -44.7365 6.08998e-05 0.000532295 1 1 46
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -48.08 -47.68 2.05 1.94 -51.83 -45.50 1.00 644 919
## b[1] 0.11 0.11 0.51 0.50 -0.71 0.95 1.01 514 928
## b[2] 1.08 1.08 0.19 0.18 0.76 1.40 1.00 2210 1318
## b[3] 1.25 1.24 0.25 0.25 0.85 1.67 1.00 1227 1494
## b[4] 1.20 1.18 0.67 0.64 0.06 2.30 1.01 498 723
## b[5] 2.34 2.34 0.85 0.85 0.94 3.74 1.01 439 577
## b[6] -1.06 -1.08 1.04 0.97 -2.69 0.74 1.01 435 642
## s 1.62 1.61 0.18 0.17 1.37 1.94 1.00 1791 1427
## y1[1] -0.03 0.01 1.67 1.69 -2.84 2.69 1.00 1792 1970
## y1[2] 1.43 1.49 1.68 1.63 -1.40 4.13 1.00 1879 1890
## y1[3] 0.35 0.34 1.65 1.63 -2.31 3.08 1.00 2066 1929
## y1[4] 5.85 5.87 1.70 1.65 3.01 8.59 1.00 2030 1971
## y1[5] -1.70 -1.69 1.70 1.73 -4.61 1.03 1.00 1948 1972
## y1[6] 3.12 3.14 1.67 1.67 0.35 5.78 1.00 2064 2012
## y1[7] 4.26 4.27 1.81 1.75 1.23 7.21 1.00 1883 1866
## y1[8] 2.77 2.77 1.72 1.68 -0.02 5.58 1.00 2046 1972
## y1[9] -0.58 -0.63 1.72 1.67 -3.38 2.26 1.00 1904 1950
## y1[10] 0.26 0.27 1.69 1.69 -2.53 3.04 1.00 1921 1486
## y1[11] 3.08 3.05 1.72 1.69 0.38 6.01 1.00 2180 1864
## y1[12] 4.84 4.87 1.72 1.72 2.01 7.61 1.00 2176 1634
## y1[13] -0.09 -0.06 1.70 1.64 -2.93 2.68 1.00 1711 1807
## y1[14] 0.81 0.83 1.77 1.74 -2.08 3.66 1.00 2223 1974
## y1[15] -2.34 -2.43 1.72 1.68 -5.11 0.49 1.00 2020 1974
## y1[16] -1.39 -1.43 1.72 1.72 -4.17 1.46 1.00 1633 1841
## y1[17] 3.87 3.87 1.71 1.65 0.94 6.72 1.00 1893 1770
## y1[18] 1.08 1.10 1.71 1.68 -1.67 3.90 1.00 2014 2015
## y1[19] 2.94 2.95 1.72 1.74 0.19 5.82 1.00 1873 1973
## y1[20] 3.53 3.53 1.71 1.69 0.77 6.32 1.00 1925 1931
## y1[21] -0.42 -0.40 1.77 1.77 -3.37 2.38 1.00 1899 1825
## y1[22] 3.40 3.42 1.80 1.79 0.39 6.37 1.00 1564 1845
## y1[23] 3.23 3.19 1.74 1.70 0.36 6.13 1.00 2135 2015
## y1[24] 2.72 2.67 1.75 1.74 -0.08 5.61 1.00 1948 1884
## y1[25] 1.10 1.12 1.73 1.75 -1.81 3.86 1.00 2053 1952
## y1[26] -1.43 -1.45 1.66 1.62 -4.29 1.30 1.00 1957 1915
## y1[27] 3.66 3.67 1.73 1.71 0.82 6.52 1.00 2069 1767
## y1[28] 1.66 1.69 1.72 1.64 -1.10 4.53 1.00 2075 1912
## y1[29] 0.90 0.92 1.80 1.82 -2.10 3.87 1.00 1575 1744
## y1[30] -2.69 -2.68 1.82 1.82 -5.63 0.25 1.00 1607 1778
## y1[31] 2.52 2.51 1.85 1.79 -0.45 5.54 1.00 1685 1845
## y1[32] -0.05 -0.08 1.77 1.75 -2.93 2.86 1.00 1968 1919
## y1[33] 0.78 0.81 1.75 1.71 -2.09 3.66 1.00 1950 1930
## y1[34] 2.30 2.32 1.74 1.67 -0.66 5.16 1.00 1937 1855
## y1[35] -1.02 -1.01 1.74 1.77 -3.88 1.85 1.00 1841 1880
## y1[36] 1.66 1.68 1.72 1.68 -1.20 4.57 1.00 1821 1390
## y1[37] 0.96 0.95 1.66 1.69 -1.75 3.65 1.00 2029 1899
## y1[38] -1.08 -1.12 1.73 1.71 -3.90 1.77 1.00 1849 1758
## y1[39] 2.27 2.29 1.70 1.72 -0.51 5.11 1.00 1918 1650
## y1[40] 2.96 3.02 1.76 1.71 0.02 5.84 1.00 1915 1933
## y1[41] -1.92 -1.93 1.71 1.68 -4.70 0.85 1.00 2231 2057
## y1[42] 3.46 3.44 1.75 1.64 0.59 6.40 1.00 1661 1934
## y1[43] 2.00 2.00 1.76 1.71 -0.84 4.80 1.00 2039 1738
## y1[44] -0.78 -0.78 1.69 1.66 -3.50 1.99 1.00 1941 1697
## y1[45] -0.34 -0.38 1.76 1.71 -3.26 2.60 1.00 1448 1773
## y1[46] 0.21 0.29 1.79 1.83 -2.78 3.06 1.00 1675 1887
## y1[47] 3.26 3.25 1.73 1.74 0.51 6.13 1.00 2209 1916
## y1[48] -3.37 -3.33 1.78 1.71 -6.31 -0.51 1.00 1566 1931
## y1[49] -1.11 -1.11 1.69 1.67 -3.80 1.67 1.00 2057 2059
## y1[50] 1.07 1.10 1.71 1.70 -1.67 3.83 1.00 1892 1837
## m1[1] -0.05 -0.03 0.53 0.51 -0.92 0.79 1.00 1094 1063
## m1[2] 1.39 1.40 0.43 0.41 0.68 2.11 1.00 1340 1244
## m1[3] 0.31 0.31 0.41 0.40 -0.35 0.98 1.00 1569 1384
## m1[4] 5.88 5.87 0.59 0.60 4.95 6.84 1.00 2809 1845
## m1[5] -1.69 -1.69 0.49 0.49 -2.49 -0.88 1.00 2168 1366
## m1[6] 3.15 3.14 0.43 0.42 2.44 3.87 1.00 2123 1609
## m1[7] 4.32 4.30 0.68 0.67 3.20 5.43 1.00 1740 1627
## m1[8] 2.74 2.73 0.57 0.56 1.78 3.67 1.00 1322 1478
## m1[9] -0.62 -0.62 0.50 0.49 -1.43 0.19 1.00 1939 1596
## m1[10] 0.29 0.27 0.52 0.51 -0.57 1.15 1.00 1756 1483
## m1[11] 3.04 3.05 0.52 0.50 2.20 3.89 1.00 1405 1458
## m1[12] 4.82 4.81 0.51 0.52 4.00 5.67 1.00 2573 1510
## m1[13] -0.07 -0.07 0.55 0.53 -1.00 0.81 1.00 1451 1544
## m1[14] 0.87 0.88 0.68 0.67 -0.30 1.95 1.00 1809 1477
## m1[15] -2.34 -2.35 0.52 0.51 -3.18 -1.50 1.00 2373 1410
## m1[16] -1.40 -1.41 0.55 0.54 -2.27 -0.53 1.01 593 1016
## m1[17] 3.87 3.89 0.59 0.57 2.92 4.86 1.00 1390 1411
## m1[18] 1.07 1.07 0.47 0.47 0.27 1.83 1.00 1516 1381
## m1[19] 2.97 2.97 0.71 0.68 1.78 4.13 1.00 1179 1262
## m1[20] 3.48 3.49 0.49 0.48 2.69 4.32 1.00 2175 1584
## m1[21] -0.41 -0.40 0.66 0.64 -1.51 0.67 1.00 1620 1430
## m1[22] 3.50 3.52 0.65 0.61 2.42 4.57 1.00 716 926
## m1[23] 3.23 3.23 0.43 0.44 2.54 3.94 1.00 2305 1679
## m1[24] 2.71 2.70 0.64 0.65 1.66 3.77 1.01 775 1206
## m1[25] 1.07 1.07 0.63 0.59 0.03 2.09 1.00 1915 1669
## m1[26] -1.41 -1.41 0.46 0.46 -2.13 -0.66 1.00 2115 1475
## m1[27] 3.69 3.70 0.59 0.57 2.72 4.68 1.00 1853 1521
## m1[28] 1.68 1.69 0.46 0.44 0.94 2.44 1.00 1260 1400
## m1[29] 0.92 0.93 0.73 0.70 -0.27 2.12 1.01 608 1039
## m1[30] -2.68 -2.69 0.69 0.67 -3.81 -1.58 1.01 507 897
## m1[31] 2.52 2.52 0.78 0.78 1.24 3.78 1.01 628 935
## m1[32] -0.03 -0.04 0.55 0.54 -0.94 0.88 1.00 1770 1569
## m1[33] 0.82 0.84 0.61 0.60 -0.22 1.82 1.00 845 1303
## m1[34] 2.27 2.28 0.54 0.53 1.39 3.12 1.00 1725 1563
## m1[35] -0.98 -0.96 0.68 0.67 -2.12 0.13 1.00 1143 1335
## m1[36] 1.74 1.75 0.64 0.64 0.67 2.78 1.00 1045 1283
## m1[37] 0.89 0.89 0.47 0.47 0.12 1.67 1.00 1690 1527
## m1[38] -1.07 -1.06 0.51 0.50 -1.92 -0.24 1.00 1808 1516
## m1[39] 2.26 2.27 0.46 0.45 1.51 3.01 1.00 1371 1414
## m1[40] 2.99 3.00 0.51 0.49 2.17 3.84 1.00 2097 1661
## m1[41] -1.91 -1.92 0.56 0.54 -2.82 -0.99 1.00 2115 1560
## m1[42] 3.48 3.49 0.57 0.54 2.56 4.41 1.00 1821 1597
## m1[43] 1.91 1.91 0.60 0.58 0.91 2.85 1.00 1645 1521
## m1[44] -0.81 -0.80 0.45 0.44 -1.53 -0.08 1.00 1869 1454
## m1[45] -0.33 -0.33 0.53 0.52 -1.18 0.55 1.01 490 946
## m1[46] 0.23 0.24 0.64 0.63 -0.86 1.28 1.00 1100 1341
## m1[47] 3.27 3.29 0.61 0.57 2.29 4.26 1.00 1690 1625
## m1[48] -3.33 -3.33 0.73 0.72 -4.53 -2.15 1.01 544 910
## m1[49] -1.07 -1.07 0.45 0.44 -1.77 -0.35 1.00 1979 1471
## m1[50] 1.01 1.01 0.41 0.40 0.34 1.68 1.00 1540 1448
fn(f2)## Initial log joint probability = -423.203
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 39 -21.4026 6.40578e-05 0.000170812 1 1 46
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -25.71 -25.33 2.16 2.04 -29.85 -22.94 1.00 690 1089
## b[1] 0.22 0.23 0.30 0.31 -0.27 0.69 1.00 650 903
## b[2] 0.97 0.97 0.12 0.12 0.77 1.17 1.00 2759 1450
## b[3] 1.06 1.06 0.16 0.16 0.80 1.33 1.01 1353 1269
## b[4] 1.59 1.59 0.40 0.40 0.93 2.24 1.00 562 713
## b[5] 1.81 1.79 0.55 0.54 0.94 2.73 1.00 421 605
## b[6] -0.96 -0.96 0.12 0.12 -1.16 -0.76 1.00 2380 1685
## b[7] -1.04 -1.04 0.66 0.63 -2.12 0.03 1.00 384 571
## s 1.03 1.02 0.12 0.11 0.87 1.24 1.00 2002 1679
## y1[1] 0.81 0.80 1.11 1.10 -1.00 2.64 1.00 2046 1918
## y1[2] 1.96 1.99 1.08 1.08 0.21 3.76 1.00 2025 1893
## y1[3] 0.77 0.74 1.04 1.01 -0.95 2.47 1.00 1918 1998
## y1[4] 3.65 3.64 1.14 1.13 1.76 5.52 1.00 1793 1860
## y1[5] -2.29 -2.29 1.10 1.06 -4.11 -0.46 1.00 2053 1860
## y1[6] 3.18 3.20 1.06 1.08 1.35 4.86 1.00 1945 1971
## y1[7] 0.77 0.76 1.23 1.23 -1.23 2.84 1.00 2126 1791
## y1[8] 1.37 1.36 1.12 1.13 -0.44 3.20 1.00 1957 2013
## y1[9] 0.09 0.09 1.10 1.13 -1.73 1.83 1.00 1968 1881
## y1[10] -0.31 -0.32 1.09 1.08 -2.09 1.48 1.00 2098 2012
## y1[11] 2.78 2.78 1.14 1.08 0.92 4.63 1.00 1931 1863
## y1[12] 3.65 3.63 1.07 1.06 1.93 5.40 1.00 1910 1895
## y1[13] -0.84 -0.81 1.09 1.03 -2.67 0.93 1.00 1695 1918
## y1[14] 3.38 3.38 1.18 1.15 1.36 5.34 1.00 1939 1893
## y1[15] -3.78 -3.75 1.11 1.10 -5.59 -2.00 1.00 2259 2020
## y1[16] -1.29 -1.33 1.09 1.05 -3.07 0.56 1.00 1871 1720
## y1[17] 2.90 2.90 1.06 1.07 1.23 4.65 1.00 1803 1751
## y1[18] 1.04 1.05 1.09 1.10 -0.72 2.78 1.00 1666 1921
## y1[19] 4.91 4.89 1.16 1.12 3.01 6.80 1.00 1933 1929
## y1[20] 3.42 3.39 1.08 1.05 1.63 5.23 1.00 1965 1933
## y1[21] 2.17 2.20 1.16 1.14 0.22 4.01 1.00 2070 1802
## y1[22] 3.13 3.10 1.14 1.13 1.25 4.97 1.00 1739 1971
## y1[23] 3.13 3.13 1.03 1.00 1.38 4.82 1.00 2011 2014
## y1[24] 1.72 1.68 1.10 1.06 -0.07 3.56 1.00 1965 1799
## y1[25] 2.12 2.08 1.10 1.07 0.38 3.96 1.00 2056 2040
## y1[26] -1.92 -1.95 1.04 1.03 -3.60 -0.18 1.00 1907 2015
## y1[27] 5.19 5.22 1.13 1.10 3.31 7.07 1.00 1921 2145
## y1[28] 2.20 2.20 1.10 1.06 0.37 3.97 1.00 1957 1869
## y1[29] 2.31 2.34 1.14 1.11 0.39 4.16 1.00 1833 1775
## y1[30] -3.49 -3.44 1.09 1.06 -5.29 -1.70 1.00 1630 1729
## y1[31] 2.64 2.65 1.17 1.24 0.76 4.48 1.00 1396 1818
## y1[32] -0.80 -0.80 1.14 1.09 -2.72 1.07 1.00 1917 1765
## y1[33] 0.38 0.41 1.09 1.08 -1.41 2.14 1.00 1699 1891
## y1[34] 3.39 3.41 1.09 1.05 1.56 5.19 1.00 2024 1910
## y1[35] -2.88 -2.86 1.18 1.14 -4.91 -0.93 1.00 1922 1921
## y1[36] 2.76 2.76 1.11 1.07 0.98 4.58 1.00 1750 1919
## y1[37] 0.60 0.58 1.05 1.01 -1.11 2.31 1.00 2019 2032
## y1[38] -0.73 -0.72 1.06 1.03 -2.50 1.00 1.00 2079 1875
## y1[39] 2.52 2.50 1.10 1.06 0.79 4.34 1.00 1962 1884
## y1[40] 3.32 3.32 1.08 1.10 1.56 5.05 1.00 1861 2034
## y1[41] -2.30 -2.33 1.10 1.08 -4.10 -0.47 1.00 1894 1893
## y1[42] 4.78 4.80 1.08 1.07 2.98 6.56 1.00 2186 1921
## y1[43] 2.33 2.33 1.09 1.05 0.54 4.15 1.00 1957 1973
## y1[44] -0.78 -0.75 1.10 1.09 -2.67 1.02 1.00 1779 1766
## y1[45] -0.20 -0.19 1.09 1.04 -2.03 1.57 1.00 2060 2010
## y1[46] 0.02 -0.02 1.12 1.10 -1.77 1.93 1.00 1961 1800
## y1[47] 5.15 5.15 1.12 1.10 3.33 6.96 1.00 1777 1856
## y1[48] -4.87 -4.88 1.12 1.12 -6.71 -3.02 1.00 1646 1777
## y1[49] -1.21 -1.25 1.09 1.07 -3.02 0.68 1.00 2111 2109
## y1[50] 1.56 1.57 1.08 1.03 -0.23 3.27 1.00 1998 1890
## m1[1] 0.83 0.83 0.35 0.33 0.25 1.39 1.00 1167 1178
## m1[2] 1.90 1.90 0.28 0.28 1.44 2.35 1.00 1107 1263
## m1[3] 0.76 0.76 0.26 0.26 0.33 1.20 1.00 1335 1310
## m1[4] 3.62 3.61 0.47 0.47 2.85 4.41 1.00 2459 1349
## m1[5] -2.26 -2.27 0.31 0.31 -2.76 -1.76 1.00 2097 1527
## m1[6] 3.17 3.18 0.27 0.27 2.74 3.62 1.00 1683 1735
## m1[7] 0.77 0.76 0.61 0.61 -0.21 1.75 1.00 2135 1509
## m1[8] 1.35 1.35 0.38 0.37 0.73 1.99 1.00 1495 1559
## m1[9] 0.10 0.10 0.32 0.32 -0.43 0.62 1.00 1673 1472
## m1[10] -0.30 -0.30 0.34 0.32 -0.87 0.25 1.00 1704 1643
## m1[11] 2.81 2.81 0.33 0.34 2.25 3.34 1.00 1282 1311
## m1[12] 3.65 3.66 0.36 0.35 3.08 4.25 1.00 2231 1573
## m1[13] -0.83 -0.83 0.36 0.35 -1.43 -0.27 1.00 1426 1481
## m1[14] 3.35 3.36 0.54 0.52 2.48 4.21 1.00 2214 1608
## m1[15] -3.75 -3.75 0.37 0.37 -4.37 -3.15 1.00 2540 1639
## m1[16] -1.30 -1.30 0.33 0.32 -1.85 -0.77 1.00 741 943
## m1[17] 2.88 2.88 0.40 0.39 2.24 3.52 1.00 1359 1430
## m1[18] 1.04 1.04 0.29 0.28 0.56 1.51 1.00 1324 1445
## m1[19] 4.91 4.90 0.50 0.49 4.09 5.70 1.00 1666 1654
## m1[20] 3.43 3.43 0.31 0.31 2.93 3.96 1.00 2078 1713
## m1[21] 2.15 2.16 0.54 0.52 1.29 3.02 1.00 1989 1588
## m1[22] 3.13 3.13 0.43 0.42 2.46 3.84 1.00 762 889
## m1[23] 3.10 3.10 0.27 0.27 2.66 3.54 1.00 1871 1815
## m1[24] 1.71 1.71 0.40 0.40 1.04 2.35 1.00 1098 1362
## m1[25] 2.15 2.14 0.43 0.40 1.47 2.87 1.00 2401 1589
## m1[26] -1.88 -1.88 0.29 0.29 -2.35 -1.41 1.00 2023 1409
## m1[27] 5.13 5.14 0.39 0.39 4.51 5.76 1.00 1703 1474
## m1[28] 2.21 2.21 0.30 0.31 1.72 2.69 1.00 985 1213
## m1[29] 2.35 2.35 0.47 0.47 1.57 3.09 1.00 813 1196
## m1[30] -3.47 -3.46 0.42 0.42 -4.17 -2.78 1.00 675 1020
## m1[31] 2.64 2.63 0.51 0.52 1.84 3.48 1.00 639 909
## m1[32] -0.79 -0.79 0.36 0.35 -1.39 -0.20 1.00 1792 1718
## m1[33] 0.38 0.38 0.40 0.40 -0.27 1.06 1.00 910 1076
## m1[34] 3.37 3.37 0.37 0.35 2.77 3.97 1.00 1586 1421
## m1[35] -2.89 -2.89 0.51 0.49 -3.73 -2.05 1.00 1254 1347
## m1[36] 2.78 2.76 0.42 0.43 2.10 3.44 1.00 1318 1653
## m1[37] 0.59 0.60 0.30 0.28 0.08 1.08 1.00 1532 1712
## m1[38] -0.74 -0.73 0.32 0.31 -1.26 -0.22 1.00 1420 1725
## m1[39] 2.48 2.48 0.30 0.30 1.99 2.96 1.00 1224 1256
## m1[40] 3.29 3.29 0.33 0.32 2.77 3.84 1.00 2119 1452
## m1[41] -2.27 -2.27 0.35 0.36 -2.84 -1.69 1.00 2052 1722
## m1[42] 4.77 4.78 0.37 0.37 4.17 5.37 1.00 1657 1451
## m1[43] 2.36 2.36 0.38 0.37 1.74 2.97 1.01 1737 1558
## m1[44] -0.73 -0.73 0.28 0.27 -1.18 -0.27 1.00 1594 1644
## m1[45] -0.19 -0.18 0.31 0.32 -0.70 0.30 1.00 616 792
## m1[46] -0.01 0.01 0.41 0.41 -0.69 0.66 1.00 1183 1323
## m1[47] 5.14 5.15 0.43 0.43 4.46 5.83 1.00 1639 1394
## m1[48] -4.85 -4.85 0.47 0.47 -5.63 -4.07 1.00 793 1255
## m1[49] -1.22 -1.21 0.28 0.28 -1.67 -0.75 1.00 1793 1489
## m1[50] 1.57 1.58 0.27 0.27 1.14 2.02 1.00 1328 1070
